



3 1176 00162 5855 


NASA Contractor Report 159342 


NASA-CR- 159342 
19800024641 


OZONE DATA AND MISSION SAMPLING ANALYSIS 


John L. Robbins 


KENTRON INTERNATIONAL, INC. 
Hampton, Virginia 23666 



OCT 3 1980 


Ai JbLc, '(-.SEARCH CENTt<» 
LIBRARV, NASA 


NASA Contract NASI -16000 
September 1980 


IWNSA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton, Virginia 23665 



NF01172 



w 


T 


OZONE DATA AND MISSION 
SAMPLING ANALYSIS 

by 

John L. Robbins 
September 1980 

V-19100/0LTR-032 







SUMMARY 


Techniques have been developed to analyze global data sets of atmospheric 
constituents and to evaluate mission sampling strategies using these global 
data sets. Mathematical formulations and computer programs were developed 
to reduce and model global data fields and to perform statistical analyses 
of results. 

The grouping scheme used to reduce data into a global grid network is shown 
and data storage methods are discussed. Procedures for modeling these data 
with spherical harmonic functions and empirical orthogonal functions (EOF) 
are detailed mathematically and numerical computer solutions are developed. 
Eigenanalysis techniques in conjunction with these EOF models are illustrated 
for reducing the dimensionality of large data sets. 

The seemingly ever-present "missing data" problem is examined using the 
sample autocorrelation function. A linear regression technique is demon- 
strated which generates a "corrected" ozone satellite data set based on 
Dobson spectrophotometer (land based) measurements. 
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I. INTRODUCTION 


Defining the temporal and spatial variability of atmospheric constituents 
requires a sampling strategy and sensing technique that is consistent with 
the nature of the species being studied. Measurements of the global ozone 
field have been made for years from the ground^, from aircraft and balloons^, 
and more recently from satellites^. This information can be examined to deter- 
mine something about the statistical nature of these data and, generally, 
the types of sampling schemes that should be considered. 

The objective of this sampling study is to evaluate various sampling 
schemes which are based on the current understanding of the global ozone 
field and on other mission related constraints. To accomplish this, repre- 
sentative data must be acquired and reduced into a usable form. A model of 
the global ozone field must be developed. Computer simulated missions can 
be generated by "measuring" the global ozone field as represented by this 
model as viewed by selected sampling schemes. How well these sampling 
missions "recover" the model is determined by statistical analysis techniques 
which serve in the mission evaluation process. 

This report addresses itself not so much to the overall sampling 
evaluation problem but to the techniques that have been thus far developed 
and utilized toward that end, especially in the areas of preliminary data 
manipulation and reduction, model development, computer simulations of 
sampling missions, and the associated statistical analysis techniques used 
throughout the work. 

Appendix A shows the primary computer programs mentioned during the 
discussion. 
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II. PRELIMINARY DATA ANALYSIS - DATA MANIPULATION AND REDUCTION 

Global ozone data utilized in this study are primarily from 
the Backscattered Ultraviolet (BUV) experiment aboard the Nimbus 4 satellite. 
These data have been received from the National Space Science Data Center in 
the form of IBM unformatted binary magnetic tapes. To some lesser extent 
ozone measurements from Dobson spectrophotometers are used. These 
data will be mentioned later in this report. This section is concerned 
with the BUV data. Particular items discussed include: 

1. Conversion of the data tapes from IBM internal format to NOS-CDC 
internal format 

2. Preliminary data analysis 

3. Data grouping 

(a) Global grid system 

(b) Statistical analysis 

(c) Data retrieval technique 
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1. IBM Format to CDC Format Conversion 


The BUV ozone data used for this study have been received on magnetic tape 
written in IBM internal format. In order to generate a NOS compatible set 
of data tapes the 32 bit IBM words must be unpacked into 60 bit CDC words, 
and the IBM internal format must be converted to CDC internal format. 

A computer program (BUVC0P2) was written to accomplish this task. This 
program has successfully generated a set of NOS tapes containing global 
ozone data from April 10, 1970 through May 6, 1977. Table B-1 shows the 
time coverages and designations of the various magnetic tapes involved 
in this process. Appendix B discusses the data tape structure and the 
IBM to NOS-CDC internal format conversion in more detail. 

2. Preliminary Data Analysis 

Once a set of usable data tapes have been acquired, they must , be carefully 
reviewed to ensure that their general format and content are consistent with 
the user's understanding and that there are no apparent problems with the 
data. A computer program (BUV3) has been written to look for particular 
problems associated with the data. These include; 

1. Out of sequence (OOS) data - data that are chronologically out of 
order. 

2. Out of range (OOR) data - a data record containing a latitude or 
longitude value out of its realistic range (-90° > latitude > 90°, 
0° > longitude > 360°), or a measured observable whose value is 
inconsistent with accompanying user information. 
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3. Inconsistent local time (ILT) data - since Nimbus 4 is in a 
Sun-synchronous orbit, the satellite should cross the equator at 
approximately the same local solar time each orbit. This is the 
case for both the ascending and descending portions of the orbit. 
However, only the ascending portion of the orbit is of concern here 
since the descending portion of the orbit is on the dark side of 
the Earth and the BUV experiment only works in the sunlight. Local 
solar times can readily be calculated by the expression, 

tg = tg + ^/15“, 

where tg is the local time, tg is the Greenwich mean solar time 
(GMT) of the observation, and is the longitude of the observation 
measured eastward from the prime meridian (PM). The analysis program 
calculates tg for observations within 5“ of the equator and compares 
them to the known local crossing time, t|^. If the difference 
te - tk is less than some predetermined acceptable At, agreement 
in the two is assumed to be good. Due to orbital considerations 
t|^ may change slightly as a function of time over several years. 

4. Repetitive data - two or more data records that occurred at either 
the same time or same position (the latter being consecutive measure- 
ments) but that differ in the values of other parameters. 

5. Duplicate data - a data record or records that exactly duplicates 
another data record, 

6. Reversed ground track - a series of data records that show the 
satellite ground track moving the wrong direction latitudinally . 

Such problems need to be identified and, where practical, eliminated. If 
known or suspected problem areas remain, one must be mindful of their poten- 
tial impact in further analyses. 
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Table 1 is representative of the information one may expect from the BUV3 
analysis program. This particular analysis table is for the third set of 
BUV data (BUV III). Items such as the number of files, number of records 
(observations), and mission duration give information useful for future 
analyses as well as confirming the data tapes' general structure and content. 
The diagnostics such as the quantity and nature of abnormal ozone values, 
help in determining what, if any, data editing must be performed. For 
example, an observation near the equator whose calculated local time, t^, 
disagrees with the known local time, tj^, significantly could mean that 
either the GMT or longitude are incorrect. However, if several observations 
near the equator for a given orbit show disagreement, the entire orbit is 
suspect and requires more careful scrutiny. Appendix C describes a linear 
approximation for calculating local time as a function of latitude that is 
used between ±60° for this purpose. 

3. Data Grouping Scheme 

The most recent set of BUV data tapes covers the period from April 10, 1970 
through May 6, 1977 and contains 1,034,456 total ozone observations. There 
are 20 parameters associated with each observation as described in Table B-2. 
This amounts to 20,689,120 computer words of data that are contained on these 
tapes. In order to work with such large quantities of data they must be 
grouped in a manageable form and stored in such a way as to be easily re- 
trievable. 

It was decided to group the data according to a global grid system each 
element of which would be 5° in latitude by 15° in longitude. This arrange- 
ment lends itself nicely to the format of a data array dimensioned 36 x 24 
where there are 36 rows representing the 5° latitudinal zones and 24 columns 
representing the 15° longitudinal sectors. This global grid system is illu- 
strated in Figure 1. The indices shown in the figure follow from the 
expressions. 
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0 ) 


(6/5°) +1, for . 0° < 6 < 90° 
(0/5°) + 19, for -90° < 6 < 0° 


and 

j = (4>yi5°) + 1, for 0° < < 360°, (2) 

where e is the latitude and <{i^ is the longitude measured westward from 
the PM. 


As the data are being grouped into this grid, it is convenient to compile 
a set of elementary statistics describing the data's behavior. Useful 
quantities that can be readily calculated for a given time period include: 


1. Sampling distribution 

2. Data means 

3. Data variances. 


The basis for this analysis is the grid format described. Each obervation 
is placed in a grid block based on its latitude (e) and longitude (cfi^) 
according to equations (1) and (2). The global sampling distribution can 
readily be determined by counting the accumulation of observations into each 
block (i,j). The zonal data distribution is found by summing this result 
over j (1 < j < 24) for each individual zone (1 < i < 36). 


For each grid block the mean ozone value, the mean position of observations, 
and the mean time of observations are calculated as shown below; 




E 

£=1 




(3) 


where the d represent the data record of either latitude, longitude, 
time, or ozone value contained in block (i,j); k.. is the number of 

I J 

observations contained in block (i,j); and X.. is the block mean for 

* V 

whichever of the above quantities is represented by d^. 
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Zonal means are calculated by 

24 


24 


(. = Z { E d ) ./ E k. .. 

^ j=l £=1 ^ ^ j=l 


( 4 ) 


or 


24 


Xi = jf, Xij 


(5) 


where 


24 

h. = E k... 

^ j=l 

Associated variance calculations follow from 


( 6 ) 


VAR(x) = <(x - <X>)^>. 


(7) 


Then the variance of the data contributing to the grid block mean becomes 


k. . 

1 j P P 

E d: - k.. xf, 

a 1 J 1 j 


XiJ •‘ij-'UM 

and the variance of the data contributing to the zonal mean becomes 
24 ^-j 


( 8 ) 


\ h.-l 


E ( E dj)j - h. Xj 

j = l £=1 ^ -J ^ ^ 


(9) 


The subscripts on show the number of ozone observations in the sample 
being considered. 


Finally, a "data mean" and variance are calculated which include all available 
data from the global grid. The data mean is 

36 24 '^ij 36 24 

X = E E (E dj../ E E k.., (10) 

i=l j=l £=1 ^ i=l j=l 
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or 


36 24 

X= E E X.. k../M. (n) 

i=l j=l 

where 

36 24 

M = E E k... (12) 

i=l j=l 

This "data mean", X, is not referred to as a global mean since the spatial 
distribution of the BUV data is non-uniform and, therefore, X is necessarily 
area biased. In addition, these data do not provide global coverage due to 
orbit and sensor design. In fact, BUV annual coverage extends only from 
approximately 80° south latitude to 80° north latitude. Otherwise, the 
extent to which the global grid is filled depends on the length of the time 
interval being considered and upon the actual portion of the BUV mission 
being examined. The latter is due to the fact that the data density per unit 
time decreases in the later years of the mission. 

The variance of the data contributing to the data mean is 



36 24 S-j 9 
E E ( E dp . . - MX 

i=l j=l £=1 ^ 


(13) 


The computer program 0ZSTAT2 was written to perform these analysis tasks. 
Graphics capabilities included in the 0ZSTAT2 program provide for each case 
a plot of the zonal means with +la^ error bars, a scatter diagram of the 

ozone distribution as a function of latitude, and histograms of the data 
sampling distribution as a function of latitude or longitude. Examples of 
the graphics output are shown in Figures 2 through 4. A listing of this 
computer program and accompanying subroutines is included in Appendix I. 


A means of storing and accessing these reduced data for specified time 
intervals is required. Typical time periods examined in this study include 
seasonal (90 days), monthly (30 days), weekly (7 days), and, less frequently, 
daily intervals. It was, therefore, decided to store this information on a 
daily basis in such a way that data for larger time intervals can conveniently 
be generated by accumulating the appropriate daily values. 
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Specific quantities that must be accessible on a daily basis per grid 
block are, 

1. Sampling Distribution 

2. Ozone Mean 

3. Average Latitude of Observations 

4. Average Longitude of Observations 

5. Average Time of Observations 

6-9. Variances Associated with Items 2-5 above. 

Rather than storing these specific nine pieces of information per grid block, 
it was decided to save the sums and the sums of the squares of the ozone, 
time, latitudinal, and longitudinal values along with the sampling distri- 
bution from which the required means and variances are readily calculable by 


= 


1J 


(14) 


and 


•^ij-1 


I -(j: X 
ij ij£ ij ij«. 


where x is the mean, is the associated variance, z x 


(15) 

is the sum, 


ij ij£ 


A/ I J I J A, 

and z x^ is the sum of the squares of the quantity for grid block (i,j). 
ij ij£ 

The £'s signify the following: 

£ = 1 Ozone, 

£ = 2 Time, 

£ = 3 Latitude, 

£ = 4 Longitude. 

The number of samples per block (i,j) is k^y 
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It was further decided to store this information on a mass storage random 
access (MSRA) file, primarily because this approach minimizes the computer 
storage problem inherent with these large data sets and also because of 
the convenience associated with utilizing the MSRA file for this kind of 
storage and retrieval process. A set of subroutines have been designed to 
access this MSRA file returning to the calling program a data array in the 
form of the standard 36 x 24 global grid system containing one of the nine 
quantities mentioned above for a given day or collection of days. These 
subroutines can be easily incorporated into computer programs requiring 
these grided ozone data without drastically affecting the program's storage 
requirement. Details concerning the MSRA file, its creation and its access 
are contained in Appendix D. 

The preliminary data analysis concepts discussed above are beneficial for 
the following reasons: 

1. Setting up a standard grid network as outlined establishes 
a basis for data analysis and lends itself nicely to making 
preliminary statistical calculations. 

2. The preliminary statistical analysis shows how the data 
distribution varies as a function of latitude and longitude 
which helps in the development of mission sampling strategies. 

3. Large data sets become more easily manageable when described 
by a global grid network which can be put into the form of a 
data array in the computer and saved on MSRA files. 
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III. STATISTICAL MODELING AND ANALYSIS TECHNIQUES 


An essential part of this mission sampling study is the development of models 
which describe the variability of the global ozone field and the statistical 
analysis techniques which can be used to evaluate these models and the sampl- 
ing schemes that they represent. The model primarily used in this work has 
been the Spherical Harmonic model, though the modeling of data with Empirical 
Orthogonal Functions has also been investigated and used to some extent 
throughout the effort. These models and certain statistical analysis 
techniques have been incorporated into computer programs which will be 
discussed. 

Cases arise where it is desirable to have a completely filled global grid 
system. The BUV data does not provide this required global coverage. A 
computer program has been prepared to handle this missing data problem using 
either a Spherical Harmonic model or a "model" based on autocorrelation 
functions. The data fill problem is discussed later in this section. 


4. Spherical Harmonic Model- Parameter Estimation and Evaluation 

The form of the spherical harmonic model chosen for this study is, 
M M 


y(e.,((>.) = z z 

* * r\ irt— I 


m=o n=m 


'■ mn mn' i ' 


D Y“ 
mn mn 




where 


Y®^(e,<}.) = cos(m(|)) 
Y°^(e,(^) = sin(m 4 >) 



for m=o 
for m>o 


(16) 


(17) 

(18) 

(19) 
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p'!J(cose) are the associated Legendre functions of degree n and order m, 

3nd ^ and are the coefficients associated with the functions Y 

mn mn mn 

respectively. is the Adolf Schmidt seminormal- 
mn ^ mn 

ization constant.^ e^. is the error associated with the i;Ui^ observation at 
colatitude e^. and longitude . 

For a given data set coefficients for a spherical harmonic model of specified 
degree and order are determined by a least squares solution that minimizes the 
sum of the squares of the residuals e^£. 

Equation (16) above can be rewritten as 

y(e,.*,)= (20) 

where both the odd and even functions, Y°^(e»<}>) ^nd Y^^(e,<})), are 
included in f( 0 ,(()), and, similarly, the coefficients and are 
included in b. Some care must be exercised in maintaining the proper ordering 
of the terms in equation (20). Note that N in equation (20) is the number 
of coefficients (and therefore the number of functions) contained in the 
model and not the degree of the model. Generally, the order and degree of 
the spherical harmonic models used in this study are equal. If a specified 
model is of order M and degree M, then 

N = (M + 1)^. (21) 

Equation (20) can be written in matrix form as 

Y = I B + E. (22) 

The double underline signifies a matrix quantity while a single underline 
denotes a vector. To minimize the sum of the squares of the residuals the 
quantity 

SS = i^E = (Y - F B)^(Y - F B) (23) 
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must be differentiated with respect to This leads to the so-called 
"normal equations" which can be solved such that 

i = fh ^ (24) 

The estimated coefficients contained in the § vector are unbiased since 
B = B . ^ 

Information regarding the sampling can also be gained from equation (24). 

To that end calculate the covariance of § as follows. Rewrite equation (24) 

B = G Y , (25) 

where G = (£^£)"^F^ is a function of sampling position only and is therefore 
treated here as a constant. The covariance matrix for B can be found by 


the law of propagation of errors^ such that 

Covar(B) = G Covar (Y) G^. (26) 

It is assumed that all components of are independent, 

Covar(y^. ,yj) = 0 for i ^ U (27-a) 

and have the same variance, 

Var(y.) = , (27-b) 

so that 

pvar(Y) = I (27-c) 

where is the identity matrix. 

Then equation (26) may be rewritten as 

Covar(B) = G i G^ (28-a) 

= [(F^£)"^ [(£^£)'^ (28-b) 

Covar(B) = [F^F)'^ F[(F^F)“^]^ o^. (28-c) 


Now consider some symmetric matrix z. Since the operations TRANSPOSE (T) 
and invserse (-1) commute^, 



But ^ IS also symmetric; therefore, 

L = ^ ^ (30) 

and 


(l"b^ = . • (31) 

As is also a symmetric matrix, by applying equation (31), equation (28-c) 
may be written as 


£ovar(B) = (e'^'i)"^ (32) 

where the variances associated with the estimated coefficients b^ are the 
corresponding diagonal elements of the covariance matrix. The off-diagonal 
elements are, of course, the covariance terms. In this study a has typically 
been set equal to one, so that 

£ovar(B) = (E^£)'\ (33) 


This is an important result that statistically describes how well the 
model can be fitted to the sample space being considered. Recall that this 
result is independent of the observation vector Y^. 


An interesting, if only heuristic, illustration is the case where only one 
sample position is contained in the sampling scheme for a spherical harmonic 
model. Consider the product of the observation matrix £ and its transpose 
written as 










P 

I 

i=l 




, (34) 


P 

Z 

i=l 


^•N^ii 
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where P is the number of observations and the function f is the same as 
it was in equation (20). Then 





(35) 


or 

t f^(6p.*p) fj(6p.*p) . 

The one sample position occurs at 

0 “0*1 “02“ . . . ~^p 

and 


( 36 ) 


*}> “ <t>1 “ 4*2 ~ *^p* 

Equation (36) then becomes 




(37) 


and equation (34) becomes 


S = F^F = P 


fl fif2 




fof 


2 ' 1 


fl 




f‘‘ 

• ••Ik 


(38) 


Now if for some matrix z all elements of a row (or column) may be obtained 
from the elements of another row (or column) by multiplication by a constant, 
that is, if z.. = (constant) z . for all j or z.. = (constant) z.. for 

I J X-J I J I X, 

all i , then det ^ = 0, ^ 


Also the inverse of a matrix can be calculated element by element accord- 
ing to 
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(39) 




CO factor ) 
det(^) 


where det(£) is the determinant of the ^ matrix. 

It can then be seen from equations (38) and (39) that 

det(£) = 0 , (40-a) 

hence the covariance matrix 

-V CO , (40-b) 

or the variance associated with estimated coefficient values would be 
infinite. 

This result demonstrates the inability of the least squares technique to 
accurately estimate the required coefficients of a model with N functions 
(N > 1) when the sample space consists of only one point. A more general 
comment that may be inferred from this example is that the variance in the 
estimated coefficient vector is a function of the sampling distribution and 
is not necessarily dependent on the number of samplings. 


5. Statistical Analysis of Spherical Harmonic Model 

2 

The data variance of the observations contained in the vector is 
calculated by 


2 _ 1 
^d ■ Tp^ 


P 

[ ^ 
i=l 


P 2 

( E yJ^/P] , 
i=l ^ 


(41) 


which comes simply from the definition of variance as in equation (7), where 
P is, as above, the number of observations. 

The RMS residual between the measurement and the spherical harmonic model is, 


RMS = 


E^E 


1/2 


(42) 


By equation (23) , 



where is the vector of estimated coefficients that minimizes the RMS. 
Equation (43) can be expanded so that 


RMS = 


1 (Y^Y - Y^F B - B^F^Y + B^F^F 
P 



(44) 


Note the term 


i = 1^0 X P) X £(P X N) X B(N X 1 ) 
is a scalar. Therefore, 


i = (I^£ i)^ = iYi 


and 


RMS = 


1 (Y^Y - 2B^F^Y + B^F 


^F B) 


1/2 


(45) 


(46) 


Substitution of equation (24) into the last term of equation (46) leads to 
RMS = (Y^Y - 2 bVy + (47) 


In this manipulation it must be remembered that 
= I, 


where is the identity matrix. 
Equation (47) quickly simplifies to 


or 


RMS = 


RMS^ = 1 (Y^Y - b'^'l^Y). 

p __ 


/2 


(48-a) 


(48-b) 


2 2 
The RMS value above is also known as the error variance, a . This is the 

e 

portion of the data variance not explained by the model. The model variance 
is then 

2 _ 2 2 

a - a , - a 
m d e 
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The ratio 

d2 - 2,2 

R - a /a . 
m d 

is often used as a criteria to judge the adequacy of the assumed model where 

0 i < 1 . 

2 

R must approach unity for the model to account for the data variability. 

The importance of terms in a given model can be measured by the relative 

power of their coefficients. Of specific interest is the power of the 

coefficients of degree n. This quantity is referred to as the degree 
2 

variance, a , and is defined as the average square of the spherical 

' Q 

harmonic solution, y^(e,(()), for degree n, or 


2tt 




2 _ 6=0 <|)=0 
°n 


(49) 


da 


9=0 4>=0 

where da is the differential area sineded(|), and 

\ " “mnC 

m=0 

The Y® (e,*) and Y° (e,*) are spherical harmonic functions 
mn mn 

as defined in equations (17) and (18). The and are the coef- 
ficients associated with Y^ (0,<)>) and Y° (6,<j)), respectively. 

rnn inii 

Letting the notation / denote the double integral '' 


(50) 


0 ,( 


, the 


6=0 (}>=0 


numerator of equation (49) may be written as 
I 


n 

m=0 6, 




(51) 
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or 


m-0 e.<j. 


mn 


e.^ 


e,(j. 


These integrals are evaluated in Appendix E so that 


4tt 


d1 <5* ]. 


^1 " 2n + 1 f « ^^nn “mn “mo 
m=0 

The denominator of equation (49) is 


Ip= I 


da = 4ir, 


e,(|) 


such that. 


0/<2" * ’)• 


or acknowledging the fact that D_ = o for m = o, 


mn 




n 2n 


W z (A^ + ). 

+ 1 ' mn mn' 


m=o 


The total power in the model coefficients can be found by summing 
M degree variances such that 


M 


M 


Total Power = z = z ^ 
n=o " n=o 2" 


W E (A^ t ) 
+ 1 ' mn mn 

m-o 


Also of interest is the integral 


>3 = 


e,(|) 


y^(e,4>) yj^(e,(!))da, 


(52) 


(53) 


(54-a) 
over the 

(54-b) 


(55) 
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or 


m=o e,4i e,i)) 


+ DA j 
mn ma„-* 


e,(|. 


Y° Y® da + D D J Y° Y°„ da] 
mn mu mn m£„''/mn m£ 

e»9 


( 56 ) 


These integrals are also evaluated in Appendix E, so that 


'3=“- 


Then the degree covariances. 


a ^ 
n£ 




<fr) da 


da 


e,4> 


are zero. 


(57) 


2 

The contribution of the zonal coefficients to is easily determined from 


Pzn = i X 100%. 

% 


(58) 


where is the percentage of the zonal contribution to the degree variance 


at degree n, and where by equation (54-a) for m = 0, 
A^ + 

7 _ on on 
n “ 2n + 1 


(59) 


is the zonal contribution for degree n. 


But does not exist for m = 0 since by equation (18 ) 
a2 

7 - ^on 

n " 2n + 1 • 


m=o,n 


( 


0,4.) = 0, so 


(60) 


20 



Substituting this result along with equation (54-a) into equation (58) gives 


on 


zn 


n 

z 

m=o 


(A, 


mn 


) 

mn' 


X 100%, 


(61) 


again remembering that = o when m = o. 

This result is useful in determining the relative importance of the zonal 
contribution to the n;^i^ degree variance. 

The model statistics discussed above, degree variance, and total power, 
explain the distribution of power in the spherical harmonic model. 

Computer program GLSRAN2 performs the various calculations mentioned thus far 
in this section. Comments concerning the associated Lengendre function recurrence 
relations utilized in GLSRAN2 are given in Appendix F. Specific details con- 
cerning file manipulations, calculation methods, and output are elaborated 
on in Appendix G. 

Further statistical analyses are performed utilizing the results of computer 
program GLSRAN2 mentioned above. These are the zonal and global means and 
variances as based on the least squares fit to the spherical harmonic model. 

First, it is desired to derive an expression for the ozone value as a 
function of colatitude only which will serve as an estimate of the zonal 
mean, 1(0). To accomplish this the model estimate as given in equation 
(16) is integrated with respect to longitude such that 

■2tt 

y(e,<f))d(}. 

i(e) = . (62) 

difi 

( f )=0 


The numerator may be written as 


■2it 

y( 0 ,(|))d<f 

( f )=0 


M 

z 

m=o 


M 

z 

n=m 


FLPn(cos0) 
mn n ' 


’2 IT 

[AmnCOs(m((,) + D^^^sin (mcj)) ]dtj> , 

( f )=0 


(63) 
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or 


‘2n 


4>=o 


Then 


M . 

y(e,(j))d(j) = 2tt z A P (cose) 
n=o 


M 


2 ( 0 ) = ^ P (cose) . 


n=o 


on n' 


(64) 


The estimated global mean is found from 


g = 


I y(e ,4>) 


da 


» 

e,<{. 


(65) 


da 


The numerator may be written as 

M M P . 

y(6,<t>) da = E ^ J P_(cose) sinede J [A cos(m<|>) 


e tp 


mn ■ n' 
m=o n-m e=o 


ii=o 


mn 


+ D^^sin(m<())] dP . 


( 66 ) 


Evaluation of the integration over (f gives 

r27T . . ( 


0 , for m ^ 0 

2iTAon» for m = 0 


(67) 


The integral over e in equation (66) may be written as 

.1 


I p[I)(cose) sinede = j P_(x) dx 
e=o " x=-i " 


( 68 ) 


using the substitution x = cose and where because of equation (67) there is 
only reason to evaluate the integral for m = 0 . 

From Appendix E 
,1 


1 P,(x) P„(x) dx = ^ 


x=-l 
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then, if«-= 0, 


x=-i 


Pq(x) P„(x) dx = 26 ^q. 


Recalling that P^{x) = 1, equation (69) may be written as 


j;, '■ 


(x) dx = 


and the numerator of equation (65) becomes 


y(e,«i>) da 


0 , 4 ' 


4ttAqo» ■fof' m = n = 0 
0 , otherwise 


Finally, since / da = 4tt, the estimated global mean is 
e,4> 





the variance of which is 
Var(g) = Var(Aj^p), 

where Var(A^jj) is calculated by equation (33). 

The global mean can also be calculated in terms of area weighted 
Another representation of the variance of the global mean can be 
from this result in terms of £ovar(£) elements. This technique 
below. 


(69) 


(70) 


(71) 


(72-a) 


(72-b) 


zonal means, 
estimated 
s developed 
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In terms of zonal means the global mean may be written us 


I 

g = s a. z./A (73) 

z i=l ^ 1 


where z^. is the estimated mean for the i^zone, the constant weighting 
factor a^ is the surface area of the i_^ zone, and 


A = z a. 
i=l ^ 


is the global surface area. Equation (73) may be rearranged as 
I 

Ag = z a. z. , (74) 

z i=l ^ ^ 


and, if the variance is taken of both sides, it becomes 

I 

Var(AgJ = Var( z a. z.) (75) 

z i=l ^ ^ 


where I is the number of zones and the z^ are to be treated as random variables. 
By the definition of variance (equation 7) the right hand side of equation 
(75) becomes 


I 

Var( z a. z.) 
i=l ^ ^ 


I 

Var( z a. z- ) 
i=l ^ ^ 


I _ I _ 2 
= <( z a. 2 . - < z a. z.» > 

i=l ’ ^ i=l ^ ^ 

= <( z z^ - E a^ <z^> )^> 
i=l ’ ^ i=l ^ 

= <( z a. l.)b , 
i=l ^ ^ 


(76) 


where 


Z^. = (z. -<zp). 


(77) 
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Equation (76) may be expanded such that 

Var( z a. 2 .) =<{a.i + (32 22 )^ + . . . 

i=l ^ ^ 

+ (aj Zj)^ + 2ai Zi 32 Z 2 + . . . + 2aj_^ Zj_ ^ Zj > 


I _ I p p M I 

Var( z a^ z.) = < I af Zt + 2 z z a. a. Z. Z,.). (78) 

i=l ^ ’ i=l ^ ^ j=l k=j+l J K j K 


Notice that 


<z^> = <(z,- -<Zi>)S= Var(z,) 


(79) 


and 


<^j ^k>= -<7j» (7k -<7k» = CovaKzj.z,,); 

then equation (78) becomes 

I _ I 2 - 

Var( z a. 2 .) = z af Var(z.) 
i=1 ^ ’ i=l ^ ’ 

I-l I 

+ 2 z z a^- a. Covar(z.- ,Zi,) . 
j=l k=j+l ^ " 

The left side of equation (75) is 


(80) 


(81) 


Var(Ag^) = Var(g^) . 


(82) 


Equating equations (81) and (82) it is found that 


I 


I-l I 


Var(g^) = 


z a. Var(z-) + 2 z z a- a^, Covar(z.,z.) 
i=l ^ ^ j=l k=j+l ^ J 


(83) 
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Var(z^.) is the variance of the mean for zone i. The colatitude position of 
zone i is taken to be at so that from equation (64) 

M 

Var(z.) = Var[ e P^Ccose.)] (84) 


or 


Var(z,) = Var( r P , A„„) (85) 

n=o 

where is the njUi^ degree Legendre function evaluated for e^- . Equation 
(85) may be evaluated in an analogous fashion to the technique used in the 
evaluation of equation (75) where is to be treated as the random 
variable. Then equation (85) may be rewritten by comparison with equation 
(83) as 

In order to complete the evaluation of equation (83) Covar( Zj ,Zj^) must be 
written in terms of known quantities. Recall that the covariance is defined 
as 


Covar(x.y) = <(x - <x» (y - <y» > . 


(87) 


Then substituting 


M 

- n^o 


( 88 ) 


into the covariance definition, equation (87), 

Covar(zjJ^) =<( P„jA„„ - P„nA„„ - P„;,A^„» 


n=o 


n=o 


n=o 


M 




M 
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Let 




( 90 ) 


then 


Covar{zj.Z|,) =<(j P„/„) ( j P„kW„) > . 


n=o 


n=o 


and 


Covar(Zj,Z|^) = <J^ '’nj''£k«„"z>- 


or 


Covar(Zj.,2,^) 


M M 

= E Z P .P„, <W W„>. 
n=o £=o ^ 


(91) 


However, by equation (90), 

<«n“z>=«*on -<^) (V -<V»>= Covar(A„„.^). (92) 

Then substituting equation (92) into equation (91) the required covariance becomes 


_ _ M M . . 

Covar(Zj.z„) = E E P Covar(A„„,A ). 
n=o £=o 


(93) 


With the help of a little algebra equation (91) may also be written as 


M 


M-1 M 


Covar(Zj.Zi^) = ^E^ ' ''Aj'’nk> <“n“z> ■ O'") 


By equation (90) 


<«n> = <(V-<\„»'>=Var(AJ 


(95) 
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Substituting equations (92) and (95) into equation (94) gives 


Covar (z Z|^) = £ Var(A^„) 

n=o 




( 96 ) 


This is a reassuring result since it reduces to equation (86) for the 
zonal variance when j = k. 

Computer programs GLOBZON and ZONVAR have been prepared to perform these 
calculations based on the results of the least squares fit to the spherical 
harmonic model, specifically the model's zonal coefficients, A , and the 
zonal elements of the covariance matrix, Covar(A^^,A^^) . When the model is 
written in the form of equation (22) these quantities are calculated from 
equations (24) and (33), respectively. 

To summarize these results the mean, z^. , for zone i as found by computer 
program GLOBZON is 

M 


where M is the degree of the model, is the n^ degree Legendre function 
for the colatitude 0^. at the center of the zone, and A^^ is the n;^ degree 
zonal coefficient. This result was shown in equation (64) and further 
developed and used in equation (84). 

The global mean, was shown by equation (72) simply to be the first 
spherical harmonic model coefficient, or 
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In order to calculate the global variance, Var(^^), the global mean was 
written out in terms of area weighted zonal means as shown in equation (73). 
The global variance was found to be 


Var(g^) = 


I 2 _ I-l I 

z af Var(z,-) + 2 z z a -a. Covar(z.,z. ) 
i=l ^ ^ j=l k=j+l ^ ^ 


as shown in equation (83). A favorable comparison of this result with 
equation (72-b), 

Var(g) = Var(A^^) , 

tends to confirm the accuracy of the zonal variance calculation as used in 
the Var(g^) calculation. The zonal variance, Var(Y^), is 


M 


M-1 M 


Var(z,) = Var(A„„) + 2 ^ 2 ^ Covar(A„„.A„^) . 


and 


Covar( Zj.z,^) 
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6. Eiqenanalysis - Empirical Orthogonal Functions 


The subject of eigenanalysis may best be introduced by means of a simple, 
if not trivial, illustration. Consider the data shown in the table below. 

Table. Data Set as Viewed in the x-| - X 2 Coordinate System. 


Observation No. 
i 

^il 

x*2 

1 

1 

1 

2 

2 

2 

3 

3 

3 


The data means and covariance matrix can be calculated as 

x-| = 2, 

52 • 2 . 


and 


COVAR 

A 




(97-a) 

(97-b) 


(97-c) 


where 


COVAR. . 

• w 






(98) 


The data variance is the sum of the diagonal terms in the covariance matrix 
or the "trace" of that matrix and is written as 


= Tr(£0VAR) = 4/3. 


(99) 
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Now with the data mean (x^ = X2 " 2) taken to be the origin of a new coordinate 
system with axes u-j and U2» the data in the table above are distributed as 
shown in the figure below. 

Figure. The u^ - U2 Coordinate System shows the "mean centered" data 

representation. 



The origin of this new coordinate system may be thought of as being displaced 
by some mean vector, m, where 

m = 2x-| + 2X2- (100) 

x-| and X2 are unit vectors along the x-j and X2 axes, respectively. 

Define another coordinate axis such that it is colinear with the data. Call 
this axis The third coordinate system is completed by placing the coordinate' 
axis through u^ = U2 = 0 and perpendicular to in the direction shown in 
Figure 5. The coordinates of the data in the i|j-] - \p2 coordinate system are 
tabulated below. 
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Table. Coordinates of Data in tjj-] - \p 2 System. 


Observation No. 

’^ii 

'f'i2 

1 


0 

2 

0 

0 

3 

/2 

0 


The covariance matrix for the data as represented in this coordinate frame is 


CpVAR = 




( 101 ) 


This analysis is of interest since it shows that the data set initially 
represented by two coordinate axes, x^ and X2» can be represented, with the 
proper translation and rotation of these axes, by only one axis, ijjp as the 
1^2 component of all three observations is zero. This result effectively cuts 
in half the amount of information required to describe this set of data. It 
follows then that the data variance must all be accounted for along the 
axis as is shown in equation (101) in accordance with equation (99). Because 
of this, and since the data mean is zero, the variance may be found from the 
mean of the sum of the squares of the i|)i axis data coordinates, or 

2 1 3 2 

a 2 . (102-a) 

i=l ^ ' 

and 



Also, by equation (101), and ip. 2 , i = 1, 2, 3, are uncorrelated since 
C0VAR(4/^ , \i/2^ ~ O' 

Now consider the case where the data set is in the form of a matrix X(MxN). 
X may be thought of as containing M measurements of an observable vector 
dimensioned by N or as N coordinate vectors dimensioned by M. The 
objective is to reduce the number of coordinate vectors required 

to accurately represent 1 and at the same time to keep account of the 


32 




data variability explained by this representation. Though more computationally 
involved, this problem is fundamentally the same as the preceding example. That 
is, by the proper selection of another coordinate system, the data may be 
arranged with respect to its coordinate axes so that the data variance is maxi- 
mized along a smaller number of its coordinate vectors and so that the various 
coordinate vectors are uncorrelated with each other, i.e.. 


C0VAR(,|/. , <|;j) = 0, for i = 1 , 2, . . ., N 

and i j. 


To this end the covariance matrix describing the data set must be diagonalized 
(all off diagonal terms are required to be zero). The covariance matrix is 


0)VAR(^ = 1 


(103) 


where the U matrix is defined such that 


“ij “ ^ij “ ’ 


(104-a) 


and 


M 


" M ^ij^ 


(104-b) 


is the data average for the j^Wi, column of X. 

Diagonalizing the covariance matrix defined by equation (103) results in a 
new covariance matrix statistically describing the data in a new coordinate 
system or "eigenspace". This covariance matrix is of the form 


A = 




(105) 
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All off-diagonal elements are zero. The diagonal elements of 4 are eigenvalues , 
or characteristic values as they are sometimes called. Associated with each 
eigenvalue is a principal axis , a coordinate axis in eigenspace. Any vector 
jjj^, as defined in equation (106) below, that is parallel to a principal axis 
is called an eigenvector . The eigenvalue equation is 

(£0VAR)j|/ = A , (106) 

which may be rewritten as 

[(£0VAR) - lA]i = 0 , (107) 

where 1 is the identity matrix. 

It is necessary to find non-trivial solutions for equation (107), that is, 
solutions where ± 0. Since equation (107) is a representation of N homo- 

geneous simultaneous equations, it can only be solved if the determinant of 
the coefficients vanishes, or 

|£0VAR - 1A| = 0 . (108) 

This is often referred to as the secular equation. Values for the scalar 
constant a which come from the solution of the secular equation are the 
sought eigenvalues. These eigenvalues are arranged in decreasing magnitude 
along the diagonal of a in equation (105). 

Once the eigenvalues are known, the associated eigenvectors can be found by 
equation (107). The N eigenvectors that pass through the origin are the 
coordinate axes in the eigenspace coordinate frame. The coordinates of 
the data in eigenspace are given by 

C = U (109) 

where ^ is defined by equations (104) and £ is a square matrix containing 
the N eigenvectors by row. The coordinates of the data in the original 
coordinate system can be found by 

1 = n + a (no-a) 

where 

^ , ( 110 -b) 
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and the matrix a, containing the N column means of X, is given by 


M ^ 


M 

z 

i=l 




(110-c) 


for k = 1 , 2 M. 

It will now be of interest to return to the earlier illustrative example solv- 
ing it from the point of view of an eigenvalue problem as developed above. 

From the data in the table (Data Set as Viewed in the x^ - ^2 Coordinate System) 
and with equations (97) and (104) the U matrix may be written as. 



( 111 ) 


This is the "mean centered" or "zero mean" data representation as shown in 
the figure above. Then by equation (103) 


CPVAR(X) = 




( 112 ) 


which is in agreement with equation(97-c). The required eigenvalues can be 
found with a little algebra and equation (108) as follows: 
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Let 



so that 


1 - X' 1 

1 1 - X' 


= 0 . 


014 ) 


(115) 


The determinant on the left hand side is readily evaluated giving 

X'^ - 2x' = 0 . (116) 

The solution of this quadratic equation is 


= 2, (117-a) 

and 

4=0. (117-b) 

or, by equation (114) , 

= 4/3, (118-a) 

and 

X2 = 0 . (118-b) 

Substituting this result into equation (107) yields 


'P 


11 


= iP 


12 


for the first eigenvalue, and 


(119) 


’^21 ^ ■’^22 

for the second eigenvalue. Here is the component of the ith eigenvector 
along the Uj axis. 
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The eigenvector associated with the first eigenvalue is any vector which has 
equal components along the u-j and U 2 axes. Then the principal axis can be 
taken as the eigenvector that passes through the origin of the u-j - U 2 coor- 
dinate system, such that the unit vector along this principal axis is 




^ ^2 

n 


( 121 ) 


where e-j and 62 are unit vectors along the u-| and U 2 axes, respectively, and 
the M'Jl is a normalization constant. 


Similarly, the unit vector along the second principal axis is 



( 122 ) 


Equations (121) and (122) can be combined and represented in matrix form as 


1 






/2 


(123) 


Here the matrix ip contains in rows the two eigenvectors that represent the 
principal axes in eigenspace. 

It can quickly be shown that i/i-i and 1 P 2 form an orthonormal set since 



From this result it follows that 


i] • ii = i2 ’ ^2 " ^ 

demonstrating that and ^ are normalized and that 
' ±2 " ^2 ' ^ 
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showing that and jj ^2 ^re orthogonal to each other. 

To digress a bit it is interesting to note that the matrix ^ in equation (123) 
is, in fact, the transpose of a rotational transformation and can be calcula- 
ted by the perhaps more conventional technique illustrated in the figure below. 

Figure. Coordinate Axes Rotation 



In the figure the primed axes, x‘ and y' , have been rotated as shown through 
an angle e. They can be represented with respect to the original axes as 


x' = X cose + y sine 

(124-a) 

and 


y' = y cose - x sine , 

(124-b) 

or written in matrix form 


(x‘ y') = (x y) /cose - sineV 

(125) 

Vsine cose/ 
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Notice that for 0 = 45° the rotational transformation in equation (125) 
becomes 



(126) 


Now, by equations (109), (111), and (123), the coordinates of the original 
data (Table. Data Set as Viewed in the x-j - X2 Coordinate System) in 
eigenspace can be calculated as 



(127) 


The original data coordinates can be reconstructed by equations (110) combined 
as 

X = a + C I . (128) 

Substituting equations (97), (123), and (127) into equation (128) gives 



(129) 
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1 V 

2 2 
3 3 


(130) 


However, as previously stated, the ^ data matrix should be retrievable by 
utilizing data along only the axis. Equation (129) for only the first 
eigenvector becomes 





X = 




(131) 


(132) 


To calculate individual 


elements of X equation (128) may be written as 


X 


ij 


+ 


N 

Z 

k=l 


^ik '^kj 


(133) 


where the i subscript on 'V' has been dropped since the column index, j, 
determines the value of now treated as a vector. 

These eigenanalysis techniques have been used to some extent in the ozone 
sampling study. Empirical orthogonal functions (EOF) have been used in the 
development of a global ozone model. These empirical orthogonal functions 
are the eigenvectors of the covariance matrix associated with a set of gridded 
ozone data. The coefficients associated with these functions are the coor- 
dinate vectors of the gridded ozone data represented in eigenspace as found 
in the C matrix defined above. 
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This EOF ozone model will be discussed below according to the following four 
development stages. 

1. Establish an appropriate data grid system. 

2. Calculate data base for model. 

3. Develop model. 

4. Test model. 

The EOF model development is based on the assumption that there will be no 
missing data blocks in the grid system. This assumption eliminates from 
consideration polar regions where there is no BUV data coverage. A further 
consideration is whether latitudinal or longitudinal variability is being 
investigated. For latitudinal variability studies data are arranged as shown 
in Figure 6-A. For longitudinal variability studies data are arranged as 
shown in Figure 6-B. The elements of the grids are found from equation (3) 
where the i and j indices are now defined as in Figures 6-A and B. 

Three data arrays constitute the minimum data base requirement for the EOF 
model. One of these arrays contains the eigenvector matrix, another contains 
the matrix of coordinate vectors in eigenspace or the coefficient matrix, and 
the last contains the N column averages of the gridded data (Figures 6-A and 
B). With one set of these three arrays the EOF model can reconstruct the 
original data grid for some specified time period. Though the EOF model is 
time independent, by supplying eigenvector, coefficient, and column average 
arrays for several time periods a model which is effectively time dependent 
can be formulated. The EOF model data base as generated during this study of 
the BUV-I data consists of one such set of arrays per week for 50 weeks. 

Data base array sets are calculated by computer program E0FA2. In this 
program the column averages are found by equation (104-b) , the eigenvectors, ^ , 
defined by equation (107) are computed by subroutine SYMQL®, and the coef- 
ficient matrix is calculated by equation (109). These data base arrays are 
saved and maintained on a MSRA file such that model data is accessible on a 
weekly basis. 
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For purposes of this discussion, it will be assumed that the source data for 
the EOF model is arranged as in Figure 6-A. The fundamental model represent- 
ation of an ozone value in grid block (i,j) is 

*ij ' * J, ‘^ik ('3'') 


as shown in equation (133). In equation (134) 

1 < n < N 

where n is the number of eigenvalues to be used by the model and is determined 
by the percentage of the total variability, P(%), to be explained or accounted 
for. The expression showing the relationship between n and P(%) is given below 
as 


P(°^) = ^ E X|, X 100% , (135) 

o k=l 

where 

2 ^ 

= E X, . (136) 

£=1 ^ 


Also, as has been demonstrated above (equations 99 and 102), the data variance 
may be written as 


= Tr(£0VAR) 


N 

E 

j=l i=l 



(137) 


The model's time dependency is incorporated by the proper selection of the a 
vector and the £ and £ matrices from the MSRA file as discussed above. 

The model development thus far makes available only the somewhat limited capa- 
bility of calculating discrete ozone values associated with grid block (i,j). 

This capability must be extended so that ozone values for specified positions 
on the Earth's surface, within the geographic boundary limitations of the 
model's data base, can be computed. This would result in a model of the form 

OZONE = X[i|-(0), C(.|)), t, P(%)]. (138) 

To this end Fourier series representations are calculated for the required eigen 
vectors and column means as a function of latitude, e, and for the required coef 
ficients as a function of longitude, <f>. Appendix H gives a brief development of 
the Fourier series representation that will be utilized below. 
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First consider approximating an eigenvector "curve" composed of discrete values. 
These values are equally spaced along a latitudinal axis and are located at 
latitude zone centers as shown on the bottom scale of Figure 7. The BUV-I data 
modeled by this technique generally has good latitudinal coverage, depending 
on the season, from the latitude zone centered at -77.5° to the zone centered 
at 77.5°. As can be seen from Figure 7 this corresponds to an eigenvector of 
32 discrete components. For the purpose of representing an eigenvector by a 
Fourier series this figure also shows certain transition scales. The "Fictitious 
Latitude Scale" simply shows the latitudinal data range where zero degrees 
corresponds to the gridded data's southern extreme. The "Fourier Scale Range" 
shows the domain of the periodical Fourier functions which will be used to 
represent the eigenvector. 

Notice that the Fourier scale range extends slightly beyond the discrete data 
scale. As far as the Fourier scale is concerned there are 33 pieces of data, 
but due to the periodic nature of the Fourier representation the functional 
value of the first discrete data point must equal the functional value of the 
last, or 

f(0°) = f(360°) . (139) 

Then over the Fourier scale range there are 32 intervals between the equally 
spaced data so that 


Fourier Scale Range 360° 

No. of Intervals over Scale 32 


11.25°/interval . 


(140) 


Both latitude scales contain 31 intervals so that the length of either latitude 
scale in terms of the Fourier scale is 

11 .25°/interval x 31 intervals = 348.75°. 

Let K-j be the conversion factor from the fictitious latitude to the Fourier 
scale such that 


K 


1 


348 . 75 ° 

155° 


2.25 . 


(141) 
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Also let 6 be the actual latitude value, 02 be the fictitious latitude value, 

6-j be the Fourier scale value, and d^ be the discrete data point number includ- 
ing any fractional part. Then, 

= k^02. (142) 

But 

02 = e + 77.5°, (143) 

so 

0^ = tc^(0 + 77.5°) . (144) 

This expression shows the relationship between the actual latitude, 0, and 
the corresponding Fourier angle, 02 * 

The relationship between dg and 02 may be written as 


^2 “ **"2^ ^0 ” ^ ’ 

(145) 

where 


<2 = 155°/31 intervals. 

(146) 

Then by equation (143) 


0 = <2^dg - 1) - 77.5°, 

(147) 

and by equation (144) 


®1 ~ *^1 *^2^^0 ” ^ ^ ’ 

(148) 


which gives the Fourier angle in terms of the discrete data point number scale. 

From the development in Appendix H the required eigenvector may be approximated 
by a Fourier series expansion 

Q 

^(e) = + E [A cos(£K ,(0 + 77.5°)) + B.sin(£K,(0 + 77.5°))] (149) 

° £=1 ^ ^ ' 

where Q = 16, since 2Q = 32 is the number of independent discrete pieces of 
data, and where equation (144) was substituted into equation (H-6) for the 
Fourier angle. 
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The procedure for finding an approximation for the mean vector "curve" is 
quite the same and leads to 

Q 

a(e) = Eq + 2 [Ej^cos(£k^{ 6 + 77.5°)) + Jj^sin(^K^(e + 77.5°))]. (150) 

1 

where only the Fourier coefficients are changed. They are found as outlined 
in equations (H-4) and (H-5). 

The Fourier representation for the coefficients is similar to that above except 
for certain scaling differences. The BUV-I gridded data ranges on the longi- 
tude scale from the block centered on 7.5° to the block centered on 352.5°. 

The Fourier angle may be written as 

<|)^ = (j. - 7.5° . 

from which by equation (H-6) 

Q 

C(<|)) = R. + L [R-Cos(Ji((}) - 7.5)) + S sin(£(<|) -7.5))] 

° ii=l ^ ^ 

where Q = 12, since 2Q = 24 is the number of discrete independent pieces of 
data, corresponding in this case to longitudinal sectors, and where the coef- 
ficients are found again by equations (H-4) and (H-5) using the already known 
discrete values of £. 

Computer program EAMODl was prepared to implement this model and to briefly 
analyze the results. To summarize the model development above as incorporated 
into the computer model consider the problem of finding an ozone value for 
some point on the Earth's surface (e‘, 4 >') at time t. Further, P'(%) of the 
data variability is to be accounted for. 

First, a set of data arrays for the appropriate time period t are accessed from 
MSRA file. Recall these three data arrays contain the eigenvector matrix, 
the coefficient matrix, £, and the column average vector, a. The number of 
eigenvectors required to achieve the specified data variability can be deter- 
mined from equation (135). 


(151) 


(152) 
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Since eigenvalues are not saved on the MSRA file, elements from the coefficient 
matrix may be used for this task. As shown above, the eigenvalue may be 
written as 

1 ^ 2 
\ " M ^ik’ 


and equation (135) becomes 


P(%) 


n 1 M 2 

k=l ^ i=l 


X 1 00% 


Equation (154) is summed iteratively over k until 
P(%) ^ P'(%) 


(154) 


(155) 


at which point n = k is taken to be the required number of eigenvectors. 

The eigenvectors are arranged by row, and the associated coefficients are 
arranged by column as stored in their respective arrays. Each of the first 
n eigenvectors are fitted according to equation (149), and each of the first 
n coefficients (column vectors) are fitted according to equation (152). 
Similarly the mean vector is represented by equation (150). The required 
ozone value can now be computed by rewritting equation (134) in terms of 
actual model parameters, instead of the global grid indices (i,j), as 

n[P'(%)] 

X[.^(e'), C(4.‘). t, P'(%)] = a^(e') + z C^(<j>’)k 4^t^e')k 

k “ 1 


where the notation n[P'(%)] signifies that the number of eigenvectors used 
is a function of the explained data variability and where the subscript t 
indicates that the arrays come from the data base for the time interval t. 

The ozone modeling technique using empirical orthogonal functions has been 
implemented and briefly analyzed by computer program EAMODl. The modeling 
aspect has been described above. As a quick evaluation of the model's 
usefulness, selected BUV-I sampling data was used to generate the required 
model data base. Then the BUV-I ozone values associated with this sampling 
were compared with the model predictions for those values. Within the 
latitudinal limits of the model, the errors ranged from 0% to 10%. 


46 



7. Data Fill Technique by Autocorrelation Functions 

The autocorrelation function typically thought of as being associated with 
time series analysis has been somewhat modified here and has been engineered 
into a data fill technique on a spatial basis. 

R(k) = E[x^ (157) 

is the definition of the autocorrelation function where E is the expectation 
operator and the set of x^ , i = 0, 1 , 2, . . N, is "zero mean" data)° The 
sample autocorrelation function 

1 N-lM-1 

Rfj(l<) - ^ x^. + ll^j (158) 


is the estimate of the autocorrection function, where k, the lag, is 
representative of time separation. 

Consider the case of global spatial distribution instead of time distribution. 
Let k represent a lag of 5° latitudinal ly and a represent a lag of 15° longi- 
tudinally. The number of samples with respect to latitude is 



180° 

5° 


= 36, 


and with respect to longitude is 


(159) 


= 560° _ 24 

'^ 2 , -| 5 ° 


(160) 


Then by analogy to equation (158) 


R (k,£) 
N(k,£) 


1 

NlMl 


24-|£l 

r 

j=l 


36-jkl 


(161) 


However, in accordance with the latitudinal index convention as shown in 
Figure 1, equation (161) is written as 


R (k,£) = 
N(k,£) 


1 

NdTiy 


24-|2l 18-|k| 36-lk| 


(162) 
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where 


N(k,lt) = N|^N, - Nj 

and Nj is the number of grid blocks containing no data. 

Now consider the data block (ij), shown in the figure below, containing no 
data. 


























i ,j 


























i-3 i-2 i-1 i i+1 i+2 i+3 


j+3 

j+2 

j+1 

j 

j-1 

j-2 

j-3 


In a sense the objective is to find a weighted mean of the 48 blocks surround- 
ing (i,j) which will serve as the "fill-in" value for the block (i,j). 

In general a weighted mean may be written as 


n n 

X = Z a- X./ E , (163) 

i=l ^ ^ i=l 

where is the weighting factor associated with x^- . Should some x^ have 
no value, indicated by x^- = 0, over the range 1 £ i £ n, then equation (163) is 
written as 


n n 

X = Z a,. X./ Z 6,-a.- 

i=l ^ ^ i=l ^ 

where 

(l, for x^- 0 

^i " (0, for x^- = 0 ■ 


(164) 


(165) 
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Finding the value for the block (i ,j) is a two-dimensional problem requiring 
summation over latitude and longitude. Let Y.. be the required weighted mean. 

I iJ 

Then by equation (164) 


Y 


ij 


3 3 

z z 

_ k=-3 £=-3 
3 3 

z z 

k=-3 £=-3 


VlIM 


''^i+k,j+£ 


*^i+k,j+£ 


(166) 


where R||^|i£] ''5 now treated as a weighting factor and from equation (162) 


R (k,£) 
N(k,£) 


1 

NiMT 


24-|^-l 18-|k| 36-|k 
z z z 
j=l i=l i=19 


^ij ^i+|k| ,j+l£|' 


(167) 


The technique briefly discussed above is currently being used as implemented 
in computer program OZFILLl on two levels, partial fill and complete fill. 
Using the partial fill technique 1/2 of the total surrounding 48 data blocks 
must contain non-zero ozone values (an ozone value of zero implies no data). 
Also previously filled blocks are not included in this count. The complete 
fill technique is used without regard to the above restrictions. 
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IV. BUV CORRECTION TECHNIQUE - DOBSON DATA 


Ozone data as measured by the Dobson spectrophotometer have been investigated 
and analyzed in conjunction with the BUV sampling analysis. These data were 
obtained from the World Ozone Data Centre in Ontario, Canada and subsequently 
have been used to adjust the BUV data as will be briefly explained below. 

For a given Dobson station certain BUV measurements are selected based on 
temporal and spatial considerations in order to calculate a linear least 
squares fit between the Dobson, y^, and the BUV, y|^, data. The great circle 
distance, s, between the Dobson station (e-j, ({> 1 ) and the BUV subsatellite 
point ( 02 , <p 2 ) is given by 

s = ^ COS 0 -I COS 02 cos|(j )2 - ij>i 1 ) , (168) 

where R = 6367.3951 kilometers is the average earth radius based on the 
Clarke spheroid of 1866. The least squares fit is of the form 

' So * 

where 3^ and 3-j are the resulting regression coefficients. 

A sufficient number of Dobson stations are utilized so that the range in 
latitudinal coverage is from approximately 75° to -45°. Both 3 ^ and 3] 
may be fit as a function of latitude, 0 , by the least squares method so that 

3o “ ciqi + ttQ 2 cos 20 (170a) 

and 

3-| = cf|i + ct ^2 cos2e . (170b) 

Then the "corrected" BUV ozone measurements, y^, as "adjusted" by the Dobson 
data may be calculated from 

y^ = Oq-j + ciq2 cos 20 + (a-|i + a^2 cos2e)yjj . (171) 
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Table 1 . 


Preliminary Analysis of the BUV-III Data 


FTN 

FILE # 

1ST DAY 

LAST 

DAY 

TOT. 

DAYS 

1ST 

REC.# 

LAST 
REC. a 

TOTAL 

REC. 


ABNORM. OZ 


lEQC 

-999. 

■ 

-77. 

OTHER 

1 

99 

126 

28 

1 

23,591 

23,591 

6 

818 

0 

0 

2 

126 

153 

28 

23,592 

47,373 

23,782 

0 

889 

0 

10 

3 

154 

182 

29 

47,374 

72,143 

24,770 

2 

1,000 

0 

1 

4 

182 

210 

29 

72,144 

97,309 

25,166 

10 

995 

0 

0 

5 

210 

238 

29 

97,310 

122,760 

25,451 

11 

993 

0 

0 

6 

238 

266 

29 

122,761 

147,467 

24,707 

0 

917 

0 

0 

7 

266 

293 

28 

147,468 

171,742 

24,275 

17 

893 

0 

0 

8 

294 

322 

29 

171 ,743 

198,572 

26,830 

86 

937 

0 

0 

9 

322 

349 

28 

198,573 

226,539 

27,967 

0 

982 

1 

0 

10 

350 

364 

15 

226,540 

240,933 

14,394 

0 

508 

2 

0 

11 

365 

392 

28 

240,934 

264,729 

23,796 

2 

745 

2 

0 

12 

393 

420 

28 

264,730 

284,198 

19,469 

1 

495 

7 

0 

13 

421 

448 

28 

284,199 

302,873 

18,675 

1 

586 

0 

0 

14 

449 

490 

42 

302,874 

326,854 

23,981 

13 

680 

0 

0 

EOF 











15 

491 

518 

28 

326,855 

345,651 

18,797 

7 

495 

1 

0 

16 

'519 

546 

28 

345,652 

367,473 

21,822 

3 

701 

0 

0 

17 

547 

575 

29 

367,474 

386,978 

19,505 

0 

623 

2 

0 

18 

575 

603 

29 

386,979 

406,530 

19,552 

2 

681 

0 

0 

19 

603 

623 

21 

406,531 

420,244 

13,714 

0 

449 

0 

0 

20 

631 

658 

28 

420,245 

441 ,153 

20,909 

4 

710 

0 

0 

21 

659 

686 

28 

441,154 

462,056 

20,903 

7 

685 

0 

0 

22 

687 

714 

28 

462,057 

483,204 

21,148 

1 

667 

1 

0 

23 

715 

729 

15 

483,205 

495,325 

12,121 , 

0 

388 

0 

0 

24 

730 

757 

28 

495,326 

520,757 

25,432 

2 

814 

3 

0 

25 

758 

785 

28 

520,758 

546,946 

26,189 

32 

871 

8 

6 

26 

786 

813 

28 

546,947 

571,915 

24,969 

8 

954 

2 

0 

27 

814 

841 

28 

571 ,916 

594,889 

22,974 

4 

671 

0 

0 

28 

842 

855 

14 

594,890 

607,974 

13,085 

0 

334 

1 

0 

EOF 











TOTALS 






607,974 

219 

20,481 

30 

17 


■ SUMMARY: 


NO. OF FTN FILES 

28 

NO. OF RECORDS 

607,974 

ABNORMAL OZONE 


OZ = -999. 

219 

OZ = - 77. 

20,481 

OTHER 

30 


RECORDS SUCH THAT ABSOLUTE LATITUDE < 5° 
BAD CROSSING TIMES (lEQC) 
OBSERVATIONS ON EQUATOR 
TOTAL DAYS ON TAPE 


NOT INCLUDED IN THE INTERVAL 
.200 < ABSOLUTE OZONE VALUE < .650 


34,698 

17 

0 

757 



LATITUDE (DEGREES) 


Figure 1. Global Grid System as Developed in Computer Program 
0ZSTAT2's Data Grouping Scheme 



LATITUDE INDICES 














TOTRL OZONE IflTN-Cfl) 


Figure 3. Example of Program OZSTAT2 Graphics Capability 

This scatter diagram shows the BUV ozone data 
distribution for June 22, 1970. 




NORMRLIZED NUMBER OF DATA POINTS 


Figure 4. Example of Program 0ZSTAT2 Graphics Capability 

This histogram shows the latitudinal sampling 
distribution of BUV ozone data for June 22, 1970 
Actual Number of Data Points = Normalized Number 
of Data Points x 103. 







Figure 6-A. EOF Model Arrangement for Latitudinal 
Variability Studies 


= 12 3 



LONGITUDE (DEG) 
(EASTWARD) 


LATITUDE (DEG) 


Figure 6-B. EOF Model Arrangement for Longitudinal 
Variability Studies 
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N = 24 
90 
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Figure 7. Transition Scales from the Lc .tude 
Scale to the Fourier Scale for 
Eigenvector Representation. 
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APPENDIX A - PRIMARY COMPUTER PROGRAMS MENTIONED THROUGHOUT MEMORANDUM 



Program Name 

Purpose 

Reference Section 

1. 

BUVCOP2 

Convert magnetic tapes from IBM 
to NOS-CDC internal format. 

IBM Format to CDC Format Conversion 
(section 1) and APPENDIX B. 

2. 

BUV3 

Preliminary data analysis program. 

Preliminary Data Analysis (section 2). 

3. 

0ZSTAT2 

Group data into global grid 
system. Perform elementary sta- 
tistical calculations. Generate 
statistical graphics describing 
global grid grouping. 

Data Grouping Scheme (section 3). 

4. 

GLSRAN2 

Regression and statistical ana- 
lysis for polynomial expansions, 
spherical harmonic functions, 
Fourier functions, and other 
specified functions. 

Spherical Harmonic Model (section 4). 

Statistical Analysis of Spherical 
Harmonic Model (section 5). 

5. 

GLOBZON 

Calculate global and zonal 
means based on spherical harmonic 
model coefficients. 

Statistical Analysis of Spherical 
Harmonic Model (section 5). 

6. 

ZONVAR 

Calculate global and zonal 
variances based on zonal 
elements of covariance matrix 
describing spherical harmonic 
coefficients. 

Statistical Analysis of Spherical 
Harmonic Model (section 5). 

7. 

E0FA2 

Eigenanalysis program calculates 
data base arrays for EOF model. 

Eigenanalysis - Empirical Orthogonal 
Functions (section 6). 

8. 

EAMODl 

EOF model and analysis. 

Eigenanalysis - Empirical Orthogonal 
Functions (section 6). 

9. 

OZFILLl 

Implements data fill techniques 
by autocorrelation and spherical 
harmonic functions. 

Data Fill Technique by Autocorrelation 
Functions (section 7). 



APPENDIX B 


To illustrate the IBM to NOS-CDC conversion process, the most recent set of 
data received will be considered. These data are contained on three IBM 
9-track magnetic tapes. The following information comes from documentation 
received with these data tapes. . 


1. Tape density 

- 1600 BPI 

2. Mode 

- Binary 

3. Parity 

- Odd 

4. Block (PRU) size 

- 8000 bytes 

5. Logical record length 

- 80 bytes 


All three tapes were generated on an IBM 360, and each tape contains 14 files. 

With the technique used, one physical record unit (PRU, 8000 bytes) or block 
of data is buffered into the central processor at a time. This is the equiva- 
lent of 2000 IBM words or 1067 CDC words. Figure B-1 illustrates what shall 
be referred to in the subsequent discussion as a sub-block, that is, 15, 
32-bit words arranged as eight packed 60-bit words. Sub-blocks are 480 
bits long since this is the smallest common multiple of 32 and 60. The 
conversion process is accomplished with one sub-block at a time. The pro- 
cedure as coded in Program BUVC0P2 is described below. 

The first block of data is buffered into an array A dimensioned by 1100. 
Unused storage locations of this array contain the value of zero. The first 
sub-block (eight words) from A is placed into the array C dimensioned by 

eight. The 15, 32-bit words in the sub-block are unpacked and arranged into 
15 right justified 60-bit words in Subroutine IBMWRDS. These 15, 60-bit 
words are stored in a temporary array B dimensioned by 15 and subsequently 
into the first 15 locations of an array D dimensioned by 2200. This process 
continues until all words in the data block have been stored in D. Sub- 
routine IBMFPC from the READIBM subroutine package can now convert these 
numbers to CDC internal format floating point numbers which are stored in 
an array E dimensioned by 2010. In general, the E array contains 2,010 
words . That is. 
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Number of words in E = 

Number of sub-blocks x Nunber of 32 bit words/sub-block, 

= 134 sub-blocks x 15 words/sub-blocks 
= 2,010 vwrds. 

The number of complete 20 word logical records in E is the integer part of 
2,010/20 or 100 records. The elements of the E array are finally written 20 
words (one logical record) at a time onto an output file which is stored on 
NOS 9-track tapes. This procedure for converting a block of data from IBM 
internal format to CDC internal format is shown schematically in Figure B-2. 

This process is repeated with the next block until the end of the tape 
is reached. A listing of Program BUVC0P2 and Subroutine IBMWRDS follows 
this appendix. 

A final comment regarding this conversion concerns the actual storage of 
data on magnetic tapes. The above technique converts one tape at a time. 

The program must therefore be run three times since three IBM tapes were 
received containing these data. The minimum amount of data contained on 
any one of these three tapes is 278,259 logical records or 5,565,180 
words. Since a standard NOS 9-track tape, 2,400 feet in length, will 
hold only 3,880,421 60-bit words, two of these tapes are required. Three 
NOS tapes were required to hold the data from the IBM tape with the most 
data. Tape designations, and associated coverage periods, are shown in 
Table B-1 . 

These NOS tapes have been prepared to be read with an unformatted binary 
READ, one logical record (20 words) at a time. These 20 words are listed 
in Table B-2. The six of these words stored per record on the condensed 
tapes, generated to minimize storage and reduce computer time, are indi- 
cated with an asterisk (*). 
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Table B-1 . Magnetic Tape Designations and Their Corresponding 

Time Coverages 



IBM Reel 

NOS TAPE 

NOS TAPE 

Time Period 

Designation 

Designation 

Designation 

April 10, 1970 

30906 

NV0738 


- May 6, 1971 


NV0739 


May 7, 1971 

34037 

NN1004 

NV0740 

- May 5, 1972 


NV0103 

NV0104 

May 6, 1972 

32701 

NV0333 


- May 7, 1977 


NV0334 




NV0335 


1 - Contains 20 words 

per logical record. 



2 - Contains 6 words 

per logical record - 

condensed tape. 





Table B-2. The Twenty Words that Constitute a 
Logical Record on the BUV Data Tapes 


1. Logical Sequence Number 

2. Orbit Number 

3. Year* 

4. Day of Year* 

5. Seconds of Day* 

6. Latitude* 

7. Longitude (westward)* 

8. Solar-Zenith Angle 

9-12. Monochromator N Values, (312.5 - 339.8)nm 
13-16. Photometer N Values, (312.5 - 339.8)nm 

17. A Channel Total Ozone Value 

18. B Channel Total Ozone Value 

19. Recommended Reflectivity 

20. Recommended Total Ozone* 

* Designates those six words maintained on condensed data tapes. 
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Figure B-1 . Sub-Block Structure.* 


Packed 60-bit 
Word Number 

32 BIT WORD ARRANGEMENT 

1 

32 bit word 

II 

28 bits 

2 

bits II 32 

bit word 

II 

24 bits 

3 

8 bi ts 1 1 

32 bit word 

II 

20 bits 
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12 bits II 

32 bit word 


II 16 bits 

5 

16 bits II 

32 bit word 

II 12 bits 

6 

20 bits 

1 1 32 bit 

word 

II 8 bits 

7 

24 bits 

II 32 

bit word 

1 1 ^ 

' ‘ bits 

8 

28 bits 

II 

32 

bi t word 


* Sub-block contains 480 bits of information. This is the equivalent of 
8 60-bit words or 15 32-bit words. 
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Figure B-2. Schematic Showing the Procedure Used 
to Convert a Block of IBM to NOS-CDC Internal Format 
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APPENDIX C - LINEAR APPROXIMATION FOR CALCULATING 
LOCAL TIME AS A FUNCTION OF LATITUDE 


A straight line approximation to the ascending portion of the local time 
variation for a Sun-synchronous orbit curve from -60° to +60° latitude was 
calculated and is shown in Figure C-1. The relationship between the local 
solar time, t^^, and the latitude, 6, for the observation was originally 
estimated to be 


_ 629.45 - 0 
£ 53.57 


(C-1) 


Since, selected BUV-III data, closely corroborated by TRACK2 computer 
program simulations, have led to what is thought to be a better estimate, 
that is 


o _ 604.54 - 0 
£ " 50.93 


(C-2) 


In any case, the error table shown below shows the maximum difference between 
equations (C-1) and (C-2) to be 0.1780 hours (10,68 minutes) where 

At = t° - tj^ . (C-3) 


Table. Error Analysis 


0 

At (hours) 

60° 

0.0619 

30° 

0.0909 

0° 

0.1200 

-30° 

0.1490 

-60° 

0.1780 
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Figure C-1 . 


Approximation to the ascending portion of the local time 
variation of the Nimbus 4 Sun-synchronous orbit curve. 


LATITUDE. 

DEGREES' 



LOCAL TIME, HOURS 


APPENDIX D - STORAGE OF GRIDDED OZONE DATA 
ON A MASS STORAGE RANDOM ACCESS FILE 


A global grid system in the form of an array dimensioned 36 x 24 has been 
selected to represent the BUV ozone data. Each of the 36 rows corresponds 
to a 5° latitudinal zone while each of the 24 columns corresponds to a 15° 
longitudinal sector. Associated with each of the 864 blocks of the global 
grid are nine values that must be saved and stored such that they will be 
readily accessible when needed. For each of these values there is a sep- 
arate array identified by the parameter ISET as shown in the table below. 

Table. Global Arrays Saved on Mass Storage Random 

Access File 


ISET 

Array Name 

Description 

1 

KK 

Sampling Distribution 

2 

SUMX 

Sum of ozone observations for each 
block 

3 

SUMXSQ 

Sum of squares of ozone observations 
for each block 

4 

SUMT 

Sum of observation times for each block 

5 

SUMTSQ 

Sum of squares of the observation times 
for each block 

6 

SUMLT 

Sum of the observed latitude for each 
block 

7 

SUMLTSQ 

Sum of squares of the observed latitude 
for each block 

8 

SUMLG 

Sum of the observed longitude for each 
block 

9 

SUMLGSQ 

Sum of squares of the observed longitude 
for each block 


It was decided that these arrays should be accessible on a daily basis for 
the 392 days beginning April 10, 1970 and ending May 6, 1971 or according 
to the time convention adopted during this study, NIMDAYS 100-491. 
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Making use of a mass storage random access (MSRA) file for this purpose is 
quite suitable. As can be seen, the actual data storage requirement here 
is 

-^Say F' — a?ray "" " 3,048.192 words. 

However, by specifying a particular array for a given day, or several days, 
the computer storage requirement is reduced to that needed for only one 
array plus an INDEX array mentioned below. 

This is illustrated in the following figure. 

Figure. Mass Storage Random Access File Arrangement 
of Global Data Arrays 



ISET 

MSRA Day 
No. 

1 

2 

3 

4 

5 

6 

7 

8 

9 

1 
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4 

5 

6 

7 

8 

9 

2 
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11 

12 

13 

14 

15 

16 

17 

18 

3 










4 
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390 










391 










392 

3520 

3521 

3522 

3523 

3524 

3525 

3526 

3527 

3528 


Each of the blocks (1-3528) shown in the figure represent a data array. Let 
"NDEX" be the number that specifies a particular array, and "IDAY" be the 
MSRA day number specification. Then 

NDEX = 9 X (IDAY - 1) + ISET . (D-1) 
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Since 


IDAY = NIMDAY - 99, (D-2) 

expression (D-1) may be written in terms of NIMDAY as 

NDEX = 9 X (NIMDAY - 100) + ISET, (D-3) 

For example, if the array SUMXSQ(ISET = 3) were required for NIMDAY 101, then 
NDEX = 12, 

and the 12j^ array would be accessed from the mass storage file. 

The INDEX array mentioned earlier must be present and must be dimensioned 
by (A + 1) where A is the total number of arrays on the MSRA. 

Listings of the subroutines GETDATl and GETDAT2 which access the BUV MSRA 
file follow this appendix. 
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APPENDIX E - ORTHONORMALITY PROPERTY OF SPHERICAL HARMONIC FUNCTIONS 


The functions i|J|^(x) for k = 1, 2, 3, . . . , are orthogonal over the interval 
(a.b) and, therefore, have the property that 
b 




t(;j(x)dx = 0, for i j. 


(E-1) 


If i = j, and if 
b 


[ij^.(x)]'^dx = 1, 


(E-2) 


then the functions are also normal , or normalized , and form an orthonormal 

set of functions over the interval (a,b). Equations (E-1) and (E-2) can be 

written as 
b 


a-* 


i|;^.{x) i|;j(x)dx = 


(E-3) 


where 6.., the Kronecker delta, has the property that 

• J 


'^ij 


0, for i 7 ^ j 


1 , for i = j 


(E-4) 


This concept can be expanded to include spherical harmonic functions over 
the surface of a unit sphere. Let y(e,i{)) be a function on the surface of a 
unit sphere, such that 


MM. - 

y(e,*) = z E 
m=0 n=m 


(E-5) 


where 


Z® (0,<l)) = cos(m<i)) Pjj(cose), 


mn 


and 


Z° (6,4i) = sin(m<|)) p|!J(cos0). 


mn 


(E-6a) 


(E-6b) 


The p||!J(cose) are associated Legendre functions. 
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It can be shown that 


'’p:(x)p:(x)dx={^-^ 
x=-1 

from which it follows that 


2n+l 


13 


f Pnlx) = m- ‘ 

x=-l 


n£ 


where Pp(x) and are associated Legendre functions for m = 0 or 

simply Legendre functions. 

Now consider the following integral equations which must be evaluated: 

'IT f2Tr 
0=0 (J)=0 








h = 


rZir 


^mn^®’*^^ „( 6 ,(}))da , 


6=0 (|)=0 
fZlT 




fir rCTT 

’3 ' J J 

6=0 (})=0 


The first may be written as 


P^(cose) P^(cose) cos(m(}>) sin(kcf))da 


where 


e ,(j) 


da = sinededcj) 

is the differential surface^area of a unit sphere and the notation 

’ rTT rCTl 

is equivalent to 

0 ,( j ) 0=0 ( j )=0 


(E-7) 

(E-8) 

(E-9) 

(E-10) 

(E-11) 

(E-12) 

(E-13) 


79 



Consider the integration over 41 . 


•Ztt 

cos(m({)) sin{k<t>)d 4 ) = 0 

( j )=0 


for m = k or m ^ k. Substituting this result into equation (E-12) 


or 


4(6.+ ) z°,(e,^)da = 0. 

6 ,(fi 


The next integral can be written as 


^2 = 


p’’’(cose) P^(cose) cos(m(j)) cos(ktt))da 
n ^ 


e,<i) 

or by (E-13) as 


^2 = 


P’’’(cose) P^(cose) sine cos(mij)) cos(k(f))d(j)de. 
n ^ 


6,4 


Again integrating over 4 yields 


/ 


2tt 


cos(ni!})) cos(k(}))d(}) = 


4>=0 


0 , for m k 
and m ^ 0 

TT, for m = k 
and m ^ 0 


and I 2 “ ^ m ji* k. Otherwise, equation (E-17) becomes 
fl 

I 2 = TT 


pjj(x) p‘j(x)dx 


x=-l 


(E-14) 
leads to 

(E-15) 

(E-16) 

(E-17) 

(E-18) 

(E-19) 
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where the substitutions x = cose and dx = -sinede have been made along with 
corresponding changes in the limits of integration. Substituting equation 
(E-7) into equation (E-19) leads to 



x=-l 
for m 0. 

If m = k = 0, the integral in equation (E-18) becomes 


■ 2 tt r2'n 

cos(m4)) cos{k(})) d(j> = d({i = 2ir, 
( j)=0 ( j)=0 


and 


(E-21) 


T - A 

^2 ~ 2n+l nJ- 


(E-22) 


for m = k = 0. Then the integral in equation (E-10) has been evaluated and 
can be written as 


e,^ 


I 4ir 


2n+l 


, for m = k = 0 
n£ 


( I otherwise 
' 2n+l (n-m) ! n£ mk’ 


(E-23) 


The final integral may be written as 

I3 = I Pjj(cose) p|^(cose) sin(m<}i) sin(k(j>)da. (E-24) 

e,4i 

By inspection, if m = 0, 13=0- 
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For m 0 the integration over <{) gives 


•2ti 0, for m k 

sin(m(j)) sin(k(}))d(() = 

*io "• for ra = k 

Equation (E-24) then becomes for m = k 


(E-25) 


fTT 


^3 = " 


p|I|{cose) P'^(cose) sinede, 




0=0 

which as before can be written as 

f1 


I3 = TT 


P^x) p'^(x)dx = 6 

2n+l (n-m)! 


x=-l 


Finally, the integral in equation (E-11) is 

" 2W(n^Ji '^nji^mk^mo’ 

0 ,(j) 


(E-26) 


(E-27) 


(E-28) 


where 6 is defined such that 
ab 


’ab 


0, for a = b 

1 , for a t' b 


(E-29) 


Now define 


Y® (0,d>) = F^ Z® (0,41), 
mn' mn mn'' 


and 


Y° ( 0 ,<^) = F^ Z° ( 0 , 4 ) , 
mn mn mn ^ 


(E-30) 

(E-31) 
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where 


mn 


2(n-m 

r 

(n+m' 

! 


1 , for m = 0 
1/2 

, for m > 0 


(E-32) 


It is necessary to evaluate the three integrals 






e , 4 > 




0><)> 


e ,< 


Y°„(e.^) Y°Je.^)da. 


(E-33) 


(E-34) 


(E-35) 


The right-hand sides of equations (E-33) through (E-35) may be written as 


I = fm 


0 ,<j) 


0,c 


mn 






0 ,({) 


Z“„(e,*) Z^,(9.*)da, 


0 ,(j) 


v“n(6.*) Yfc\(e.»)da= | Z^(e.t) z“^(e.+)da, 

0,({i 0,4) 

Substituting equation (E-15) into equation (E-36) yields 


(E-36) 

(E-37) 

(E-38) 


Y|;,da = 0. 


(E-39) 


0,4) 
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Similarly by equations (E-23) and (E-37) 


J Y®^(e,4,)da 2SI1 "^nAk* 

e,(j) 

Finally, equation (E-38) may be evaluated by equation (E-28) as 

(E-41) 

o »<p 

The results required for arriving at equation (53) can be found from equations 
(E-39) through (E-41), respectively. That is, 


j 

6,4' 


6,4 

and 

6,4> 



(E-42) 

[Y^„(e,*)]2da = . 

(E-43a) 

[Y:„(e.4)]'da = 

(E-43b) 


Though incidental to this discussion, it should be noted that the functions 
Y® ( 6 , 4 ) and Y° ( 6 , 4 ) are orthogonal over the unit sphere since 




= 0 


6,4 


and 


Y°„(e.*) Y“,(e.t)da = 0 


6,4 


for n ^ t, m k, or both, and 

1 O"’*) = 0 


6,4 

in any case 


(E-44a) 


(E-44b) 


(E-44c) 
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However, these functions are not normalized over the sphere as can be seen 
by equations (E-43) but are said to be semi -normalized according to Adolf 
Schmidt**by the constant defined in equation (E-32), 
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APPENDIX F - RECURRENCE RELATIONS FOR 
ASSOCIATED LEGENDRE POLYNOMIALS 


In the modeling of atmospheric constituents with spherical harmonic 
functions it is useful to have the capability of calculating the required 
associated Legendre functions using recurrence relations. Many such 
relations exist for associated Legendre polyonomial functions. 

Subroutine LEGNDR4 has been written to calculate the associated Legendre 
functions up to and including those of some specified order, m, and degree, 
n, for a given colatitude, 0. This subroutine is listed in Appendix G 
with the GLSRAN2 program. 

If P|!J(x) is the associated Legendre function of order m and degree n, then 
the first two functions are defined as^, 

Po(x)=l, (F-1) 


and 


Pf(x) = X, 


(F-2) 


where 


X = cos( 0 ) . (F-3) 

The functions of higher order and degree are evaluated by two recurrence 
relations. Consider the recurrence relation^ 

* ’) X Ox* - Ci<xn. (M) 

This expression is used to calculate zero order (m=0) functions of degree 
n+1 from the two preceding zero order terms. Setting m=0, equation (F-4) 
becomes , 

Clfx) = [(2''+’) X p“(x) - np“_,(x)], (F-5) 
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Equation (F-5) is the first recurrence relation used in subroutine LEGNDR4. 
The second recurrence relation used in LEGNDR4 conies from^ 

(2n+l){l-x^)^/2 Pjj(x) = Pj;:j](x) - P[j!](x). (F-6) 

Replacing n+1 with n and m+1 with m equation (F-6) may be rewritten as 

Pj[(x) = Pjj_ 2 (x) + (2n-l)(l-x^)^/2 

Consider the first term on the right-hand-side of equation (F-7). Since 
the order must be equal to or less than the degree of the function (see 
equation (16)), 

m £ n-2, (F-8a) 

or 

n >_m+2, (F-8b) 


and the required recurrence relation for the higher order (m>0) associated 
Legendre functions becomes. 


Pj|(x) = PQ + (2n-l)(l-x^)^/^ Pnll(^) 


(F-9) 


where 


PQ = 


p|!J_ 2 (x) , for n ^ m+2 
0, otherwise 


(F-10) 


The numerical technique described above as utilized in subroutine LEGNDR4 
has been verified up through m = n = 12 on the NOS-CDC computer system 
at NASA/LaRC. 
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APPENDIX G - THE GLSRAN2 PROGRAM 


The primary purposes of the GLSRAN2 program as used in the ozone sampling 
study are to generate global stratospheric ozone models in terms of surface 
spherical harmonic functions by performing least squares fits to sets of 
BUV data and to perform certain statistical analyses as have been outlined 
in this report ("Spherical Harmonic Model" and "Statistical Analysis 
of Spherical Harmonic Model"). The spherical harmonic model representation 
as shown in equation (20) is used by GLSRAN2. The table below shows the 
relationship between the functions, f^. , as used in this representation and 
those, Y® and Y° , as shown in equation (16). 


Table. Relationship Between Spherical 
Harmonic Function Representations 


Zonal Functions 

Sectoral Functions 

Tesseral Functions 

O O 

Q. 

II 

Ll. 

F = Y® 

V2 ni 

0 

''(3M+1)+1 " ^12 

F = p® 

F = Y° 

•^M+S ’ll 

F - 

'^(3M+l)+2 “ ’12 

F = p® 

^3 ^2 

F = Y® 

^M+4 22 

F - Y® 

‘^(3M+l)+3 ’13 

• 

F = Y° 

•^M+S ’22 

F - V® 

'^(3M+l)+4 ’13 

F = P° 

^M+1 M 

F + Y® 
■^3M ’mm 

F = Y^ 

'^(3M+1)+NT-1 M-1,M 


F = Y*^ 

'^3M+1 ’mm 

F - Y® 

'^(3M+1)+NT ’M-1,M 


In the table the functions P°, for i = 0, 1, . . M, are the zonal associated 
Legendre functions, or simply Legendre functions. M is the order and degree 
of the model . 

NT = M(M - 1) (G-1) 
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is the number of tessera! functions. There are M + 1 zonal functions and 2M 
sectoral functions. The number of terms in a model of order and degree M is 

N = (M + 1)^. (G-2) 

Model coefficients are computed according to equation (24) which may be written 
as 

B = S'^ R (G-3) 

where the "information" matrix, is defined by equation (34) and 

R = Y . (G-4) 

The S matrix, dimensioned N x N, is strickly a function of the sampling. 

As S is a symmetric matrix only its upper full triangle--the diagonal elements 
and those above the diagonal--is used in GLSRAN2. This implementation reduces 
computer time as well as the storage requirement. Since solving for the model 
coefficients requires that the inverse of S be computed, these time and stor- 
age savings become even more noteworthy. 

The upper full triangle of S is "packed" into a vector. This vector, called 1/ 
to avoid confusion, contains 

e = I (N + 1) (G-5) 

elements. The correspondence between S matrix elements and V vector elements 
is given by 

V(i) = S(m,n) (G-6a) 

where 

i = m + n(n - 1 )/2. (G-6b) 

The GLSRAN2 program is set up to either calculate the V vector based on input 
sample data or to access a previously calculated V vector through a local file. 
This is also the case for the R vector though to calculate R^ actual ozone obser- 
vations must also be available. 
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Once these data are contained on working local files, GLSRAN2 makes available 
several options regarding which model coefficients or set of coefficients can 
be computed. The S matrix elements contained on local file are associated 
with a "master" model. The most obvious option is to compute the N coeffi- 
cients for this master model. Three other options exist as listed below. 

1. Coefficients may be calculated for a model of order L (L < M). 

To do this the program selects the required "subset" of the packed 
S matrix elements contained on local file and forms a new set of 
packed S matrix elements. The same is done for the R vector. 

2. Model coefficients may be calculated based on a specified number of 
independent sampling observations (for example, a certain number 

of Dobson stations). When this option is selected the program 
determines the size of the model such that the number of model 
terms is equal to or less than the number of independent obser- 
vations and then proceeds to find the S matrix elements required 
to form the new S matrix for the subset model. 

3. Particular model coefficients may be specified according to degree, 
n, order, m, and whether they are to be associated with an odd 

(i = 1), Y°^(e,(})), or even (i = 0), Y^^(e,i}>), spherical harmonic 
function (see equations 17 and 18). Identification of required 
coefficients by this option follows from the expression: 

n + 1 , for m = 0, 

(M + 1 ) + 2m - 1 + 1 , for m = n, 

3M + m(2n - m - 1) + i, for m ^ 0 (G -7) 

and m n. 

This technique is illustrated below since the idea is fundamental to the three 
options as used to determine spherical harmonic function indices or the master 
S matrix elements required to form the subset S matrix. Assume the S matrix 
is associated with a master model of degree and order M = 5 and that the coef- 
ficients specified in the table below are sought. 
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Table. Functional Form Indices with Corresponding 
F|^ Functional Form Indices 



m 

n 

i 

k 

1 

0 

0 

- 

1 

2 

0 

3 

- 

4 

3 

2 

2 

0 

9 

4 

1 

2 

1 

18 


From the 
,e 


and 


00 

e 

o3 

e 

22 


table it can be seen that for a degree spherical harmonic model 



(G-8) 

Fg , 



Also in terms of master S matrix elements the subset S matrix for this example 
is 


'll 

=14 

^19 

^1,18 

'41 

CO 

S 49 

^4,18 

91 

S 94 

S 99 

^9,18 

18,1 

^18,4 

^18,9 

^18.18 


The following discussion pertains to the input/output (I/O) requirements and 
capabilities of 6LSRAN2. As a complete listing of 6LSRAN2 and its subroutines 
is included in this appendix the discussion is limited to I/O items involving 
the spherical harmonic model. 
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Four NAMELIST input lists control the program's operation. These are named 
below along with their associated parameters. 


1 . DATA 

(a) NDATA 

(b) MORD 

2. JOB 

(a) IDATA = 1 

= 2 

(b) IFUNC = 2 

(c) lOPT = 0 

= 1 
= 2 

(d) JOPT 

(e) ITAPE = 1 

(f) ICASE 

(g) JCASE 

3. PARAMTR 
(a) BETA 

4. J0B2 

(a) METHOD = 


(b) NFUNC 

(c) MMORD 

(d) ICODE = 


- number of observations in data set. 

- order of master model . 


- simulate a data set based on an input sampling 
scheme and model coefficients. 

- data set is an input quantity. 

- spherical harmonic model fit to be performed. 

- do not calculate S matrix. S matrix is already 
on local file TAPE4. 

- calculate S matrix and store it on local file TAPE4. 

- calculate S matrix, store it on local file TAPE4, 
and STOP program execution. 

- same description as lOPT above except that JOPT 
pertains to the R vector. 

- number of cases to be run requiring a new data set. 

- number of "sub-model " cases to be run per data set. 

- input coefficients used for data simulation. 

1 - calculate coefficients for specified subset model. 

2 - determine number of coefficients to calculate based 

on specified number of independent observations. 

3 - particular coefficients to be calculated are specified. 

4 - calculate coeffcients for complete master model. 

- number of coefficients in subset model. 

- order of subset model. 

0 - do not compute coefficients. 

1 - compute coefficients. 
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GLSRAN2 uses the FORTRAN variable dimensions source statement preprocessor 
program PRE. Variables input by this program control the size of GLSRAN2 arrays. 
These variables are: 

1. N - the number of coefficients in the master model. 

2. NN - the maximum value of NFUNC for a given run such that NN < N. 

3. NV - the number of element in the packed S matrix array such that 

NV = N(N + l)/2. 

Local files used by GLSRAN2 include: 

1. TAPEl - used for input data that must be rearranged by a user supplied 

subroutine to meet TAPE2 input file requirements. 

2. TAPE2 - standard format input data file read by subroutine REALDAT. 

3. TAPE3 - used to store such items as model coefficients and covariance 

matrix elements for future use. 

4. TAPE4 - contains elements of packed master S matrix. 

5. TAPE? - contains master R vector. 
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PROGRAM GLSRAN2 


7A Hh 


DPT-1 


FTN A. 7+485 


80/01/29. 14,11.35 


1 


5 


10 


15 


20 

UD 


25 


30 


35 


40 


PROGRAM GLSRAN2 ( I NPU T> OUTPU T> TAP E5-I NPUT# TAP E 6-QUT PUT^ 
1TAPE1*TAPE2,TAPE3,TAPE4,TAPE7) 

COMMON PI,X1#X2# YVAP# DX#X» ER # I FUNC# Y, 2 1# Z 2 # 2# M PRD» NO A T A 
DIMENS ION F( 169)#S ( 14 36 5 ) , R ( 169 ) # B ( 169 ) 

DIMENSION INDEX (169) 

DIMENSION BETA(169) 

DIMENSION NAME(15) 

NAMELIST/0ATA/X1#X2# Yl# Y2# YVAR/NDATA>21#Z2>M0RD 
NAMELIST/J0B/IDATA>1FUNC#I0PT, JOPT#ITAPE>ICASE>JCASE 
NAMELIST/PAPAMTR/PETA 

NAMELIST /JOB 2/METHOD, NFUNC/ MMORD/ ICODE 

DATA NAME/7HS »7HR ,7HC0VAR ,7HB #7HC0V 

XTHCOR ,7HZBAR »7HSY »7HRY , 7HA »7HB1 

X7HC0RI #7HRC0VAR ^7HSI »7HC0VARI / 


t 

f 


COMMENT — 

C N 

C 

C 

C 

C 

C NFUNC 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C NN 

c 

c 

C NV 

c 

c 

c 


GLSRAN2 - PARAMETER LIST. 

IS THE TOTAL NUMBER OF COEFFICIENTS ( AND THEREFORE THE NUMBER 
OF FUNCTIONS) THAT MAKE UP THE "MASTER'* MODEL. 

N IS A VAPDIM INPUT PARAMETER. 

(IF A NEW MASTER MODEL IS NOT BEING CALCULATED# N MAY BE SET 
EQUAL TO NN» SEE BELOW) . 

IS THE TOTAL NUMBER OF C QE F F 1C I ENTS ( AND THEREFORE THE NUMBER 
OF FUNCTIONS) THAT MAKE UP THE 'SUBSET* MODEL FOR A 
PARTICULAR CASE. 

NFUNC IS DETERMINED AS FOLLOWS# 

MFTHaD«l# ORDER AND DEGREE# MMQRD AND NNDEG# RESPECTIVELY# 
ARE BOTH KNOWN FOR THE DESIRED "SUBSET" MODEL. 
THEN, 


NFUNC- 1+NNDEG + MHQRD(MM0RD + 1 ) . 

METHOD-2# NFUNC-NUHBER OF INDEPENDENT OBSERVATIONS 
TO BE MODELED. 

METHOD-3# NFUNC-NUMRER OF COEFFICIENTS SPECIFIED 
TO BE MODELED. 

METHOO-4# USE ENTIRE "MASTER" MODEL. 

NFUNC-N. 

IS THE MAXIMUM VALUE OF NFUNC DURING A GIVEN RUN, 

BUT NOT TO BE LARGER THAN N. 

NN IS A VAROIM INPUT PARAMETER. 

IS THE NUMBER OF ELEMENTS IN THE V-VECTOR ( PACKED FORM OF 
THE UPPER FULL TRIANGLE OF THE S-MATRIX), 

NV IS A VaRDIM INPUT PARAMETER. 

NV-N(N+1 ) /2. 


GLSRAN2 

GLSRAN2 

PREPRQCS 

PKEPROCS 

PREPROCS 


GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 



PRHGPAM GLSRAN2 


74/7A DPT*1 


FTN A.7+A85 


80/01/29. 14,11.35 


45 


50 


55 


60 


<o 

CJ1 


65 


70 


75 


80 


C 

r 

C 

C 

r 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

w 

c 

c 

r 


NVV IS THE NUMBER OF ELEMENTS IN THE VV-VECTOP ( P AC K ED FORM OF 
THE UPPER FULL TRIANGLE OF SS-MARTRIX# S-MATRIX FOR THE 
"SUBSET” MOOED. 

NVV-NFUNC (N'=UNC + 1) / 2. 

NOTE — IN GLSPAN2, BOTH VECTORS V AND VV WILL USE THE STORAGE 
S»ACE IN THE S ARRAYCONE AT A TIME). THEREFORE# 

THE S-ARRAY MUST BE DIMENSIONED BY NV(THE LARGER OF NV AND 
NVV) WHEN THE V VECTOR IS TO BE CALCULATED. 


CMMENT - 
TAPEl 


TAPES — 


TAPE4 — 


TAPE? — 


GLSPAN2 - LOCAL FILE REQUIREMENTS. 

- DATA FOP SUBROUTINE INPUT WHICH IS TO BE REARRANGED 
AND PUT ONTO TAPE2. 

TAPE2 — LOCAL FILE CONTAINING DATA TO BE READ IN BY 
SUBROUTINE REALDAT. 

RESERVED POP RESULTS SUCH AS CALCULATED MODEL COEFFICIENTS 
HAY BE SAVED ON PERMANENT FILE OR 
PE SUBSEQUENT TO PROGRAM EXECUTION. 

CONTAIN UPPER FULL TRIANGLE OF S-MATRlX 
D THERE IN "PACKED” FORM. THIS DATA MAY 
OR MAY BE CALCULATED IN THE PROGRAM, 

CONTAIN R-MATRIX ASSOCIATED WITH THE 

F defined by the s-matrix being used. 


SO that 

THEY H 

ON MA 

GN 

ETIC TA 

LOCAL 

F 

ILE TO 

which 

I 

S STORE 

ALREA 

DY 

EXIST 

LOCAL 

f 

ILE TO 

S AMPL 

IN 

G SCHEM 


SAME 


N 

NN 

NV 

READ 

READ 


■169 
■169 
■ 14365 
(5# DATA ) 
( 5» JOB) 


GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

PREPRQCS 

PREPROCS 

PREPRQCS 

GLSRAN2 

GLSRAN2 


RFAD (5,PARAMTR) 

GLSRAN2 

WRITE (6#DATA) 

GLSRAN2 

WRITE (6#J0B) 

GLSRAN2 

WRITE (6#PARAMTR) 

GLSRAN2 


GLSRAN2 

NT — SET SEED FOR RANF AS PI*ARC COS (-1) 

GLSRAN2 

PI^ACDS (-1.) 

GL5RAN2 

CALL PAMSET (PI) 

GLSRAN2 


GLSRAN2 

DO 999 LLL-1#ICASE 

GLSRAN2 

PRINT 175, LLL 

GLSRAN2 


GLSRAN2 

IF (ITAPE.EO.l) READ(5,401) NDATA 

DATM0D2 



PRIGRAM GL5:RAN2 


7A/7A DPT-1 


FTN 4. 7+^85 


80/01/29. 14.11.35 


85 


401 

FORMAT (15) 

DATM0D2 




CALL SECOND(TIHE) 

GLSRAN2 




PRINT 300# NOATA#TIME 

GLSRAN2 



300 

FORMAT ( 1X, + NDATA-*^ I8,5X,*TIME- + #F10.3 ) 

GLSRAN2 


COMMENT — IF IQPT-0> NO NEW V-APRAY IS REQUIRED. 

GLSRAN2 

90 

C 


IF J0PT-0» NO NEW R-ARPAY IS REQUIRED. 

6LSRAN2 




IF ( JOPT.EQ.O.AND.IOPT.EQ.O) GO TO 58 

GLSRAN2 


C 



GLSRAN2 


COMMENT — INITIALIZE INPUT PARAMETERS TO GLSCDRl. 

GLSRAN2 




WbI . 

GLSRAN2 

95 

C 


F(l)-1. 

GLSRAN2 

GLSRAN2 




CALL GLSCORl ( F » S» R j W/ B# Y, N# NV# JOPT» SUM Y# YSQSUM, I ERR ) 

GLSRAN2 


C 



GLSRAN2 


C 



GLSRAN2 

100 

C 



GLSRAN2 


COMMENT — PROGRAM CHGSES EITHER TO USE REAL DATA OR TO SIMULATE 

GLSRAN2 


C 


ITSOWN DATA SUCH THAT# 

GLSRAN2 


C 


IDATA-1 DATA IS SIMULATED 

GLSRAN2 


C 


IDATA-2 REAL DATA IS READ IN 

GLSRAN2 

105 

C 



GLSRAN2 


C 


SUBSEQUENTLY THE REQUIRED MODEL FUNCTIONS ARE CALCULATED. 

GLSRAN2 




ICnUNT-0 

GLSRAN2 

VO 

CTl 


25 

CONTINUE 

GLSRAN2 




IF (IDATA.EO.l) CALL SIMDATl (BETA#N#F) 

GLSRAN2 

110 



IF (IDATA.E0.2) CALL REALDAT (F,N) 

GLSRAN2 




IF ( IFUNC.EQ.999) ST0P2 

GLSRAN2 




IF (IFUNC.E0.998) STOPS 

GLSRAN2 


C 


CALL GLSSUMl ( F # S # P # W# B > Y» N, NV# J DPT# SUM Y, YS OS U M, I E P R ) 

GLSRAN2 

GLSRAN2 

115 



lCmiNT«ICOUNT + l 

GLSRAN2 




IF (ICOUNT.EQ.NDATA ) GO TO 50 

GLSRAN2 




GO TO 25 

GLSRAN2 


c 



GLSRAN2 


c 



GLSRAN2 

120 


50 

CONTINUE 

GLSRAN2 




IF (JOPT.EQ.O) GO TD 54 

GLSRAN2 




REWIND 7 

GLSRAN2 




DO 52 I»1#N 

GLSRAN2 




WRITE (7) R(I) 

GLSRAN2 

125 


52 

CONTINUE 

GLSRAN2 




IF (J0PT.E0.2) ST0P5 

GLSRAN2 



PROGRAM GLSRAN? 


7A/7A DPT*1 


FTN A. 7+405 


80/01/29. 14.11.35 


130 


54 continue 

IF (IQPT.EO.O) 
REWIND 4 

D0_56 I«1#NV 

WRITE (4) S(I) 
56 CONTINUE 

IF (I0PT.E0.2) 
58 CONTINUE 


GO TO 58 


STOPS 


135 


140 


145 


VO 


c 

c . 

COMMENT — 
C 

c 

c 

c 

c 

c 

c 

c 

c 

COMMENT — 


FOLLOWING 
STORED ON 


STATEMENT 58# THE 
LOCAL FILES TAPE7 


REQUIRED R 
AND TAPE4# 


AND S ARRAYS ARE 
RESPECTIVELY. 


NOW USING THESE ARRAYS DR SPECIFIED SUBSETS OF THEM 
(SPECIFIED ACCORDING TO THE VALUE OF "METHOD" -SEE TABLE 
BELOW)# CALCULATE THE COEFFICIENTS OF THE REQUIRED FUNCTION 


(ACCORDING TO "IFUNC") AND STORE RESULTS 
THE JCASE PARAMETER SPECIFIES THE NUMBER 
RUN WITH CURRENT LOCAL FILE DATA. 


IN THE ARRAY "B". 
OF CASES TO BE 


GLSRAN2 

GLSRAN2 

GLSRAN2 

_GLSRAN2 

6LSRAN2 

6LSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 



C 

. METHOD 

DESCRIPTION 

NFUNC 

MMORD 

GLSRAN2 


C 



♦ ♦ ♦ ♦♦ 

♦ ♦♦♦ 

GLSRAN2 

150 

C 





GLSRAN2 


C 

1 

SPECIFY SUBSET MODEL 

NUMBER OF 

ORDER 

GLSRAN2 


C 



SUBSET MODEL 


GLSRAN2 


C 



FUNCTIONS 


GLSRAN2 


. C 





GLSRAN2 

155 

C 

2 

SPECIFY NUMBER OF 

NUMBER OF 

N/A 

6LSRAN2 


C 


INDEPENDENT OBSERVATIONS 

OBSERVATIONS 


GLSRAN2 


C 





GLSRAN2 


C 

3 

SPECIFY PARTICULAR 

NUMBER OF 

N/A 

GLSRAN2 


C 


COEFFICIENTS 

COEFFICIENTS 


GLSRAN2 

160 

c 





GLSRAN2 


c 

4 

USE COMPLETE MASTER 



“ GLSRAN2 


c 


MODEL DATA AS REQUIRED 



GLSRAN2 


G 





GLSRAN2 


C 

NOTE — 

NN IS VARIABLY DIMENSIONED 

BY THE 


GLSRAN2 

165 

C 


PREPROCESSOR PROGRAM# PRE# 

AS THE 


GLSRAN2 


r 


LARGEST NUMBER OF SUBSET FUNCTIONS 


GLSRAN2 


c 


REQUIRED IN A GIVEN RUN. 



GLSRAN2 



PROGRAM GLSRAN2 


7A/7^ OPT«l 


FTN ^.7+A85 


80/01/29. IA.11,35 


170 


175 


180 


185 


190 


to 

CO 


195 


200 


205 


DO 950 JJJ-1#JCASE 
READ (5jJ0B2) 

WRITE {6jJQB2) 

IF_ (METHOD. EQ. ^ ) GO TO 60 

IF (METHOD. EG. 1) CALL SU0S1 ( INOE X# NFUNC# MMOR D» N> MORD ) 

IF (METHOD. EG. 2) CALL SUBS2 ( INDEX, NFUNC# MORD ) 

IF (METHOD. EG. 3 ) CALL SUBS 3 ( INDEX, NFUNC, NDEG, MTiRD ) 

IF (INOEX(l) .EG.-999) GO TO 60 
GO TO 65 
60 CONTINUE 
NFUNC-N 

C NFUNC IS SET EGUAL TO N HERE FOR THE CASE OF USING THE FULL 

C MASTER MODEL WHEN METH0D»2(IE. INDEX ( 1 ) »-999 WAS RETURNED 

C FROM SUBROUTINE SUBS2). 

C THEREFORE NFUNC WILL NOT HAVE TO BE DEFINED FOR CASES 

C WHERE METHOD-A. 

DO 62 I-1,N 
INDEX(I)»I 
62 CONTINUE 

65 CONTINUE 
REWIND A 
REWIND 7 
KCOUNT-0 
JCOUNT-0 

DO 70 II-l, NFUNC 
I«INDEX( II ) 

COMMENT — IF IC0DE»0(ACC0RDING TO •J0B2* NAMELIST INPUT), DO NOT 
C COMPUTE COEFFICIENTS. THEREFORE, TAPE7 IS NOT REOUIRED. 

IF (ICQOE.EO.O) GO TO 67 

66 CONTINUE 

IF (KCOUNT.GT.INDEX(NFUNC) ) GO TO 90 

KCOUNT-KCQUNT+1 

READ (7) PR 

IF (KCOUNT.NE.I ) GO TO 66 
R(II)-RR 

67 CONTINUE 

DO 70 JJ-1,II 

COMMENT — IVV IS THE INDEX FOP THE VECTOR VV, TO BE STORED IN THE 
C S-ARRAY. 

IVV-(II^(II-l))/2+JJ 
J«INDEX( JJ) 

COMMENT — IV IS THE INDEX FOR THE VECTOR V, NOW CONTAINED ON TAPEA. 


GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSKAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 

GLSRAN2 



210 



PROGRAM GLSRAN2 


7^/74 QPT-1 


PIN 4.7+<t85 


80/01/29. 14.11.35 


ID 
ID • 





IV-(I*{I“l))/2+J 


GLSRAN2 



68 

CONTINUE 


GLSRAN2 




IF (JCnUNT.GT.IV) GO TO 90 


GLSRAN2 




JCOUNT-JCOUNT+1 


GLSRAN2 

?ri5 



READ(4) V 


GLSRAN2 




IF ( JCOUNT.NE.IV) GO TO 68 


GLSRAN2 




S (IVV) -V 


GLSRAN2 



70 

CONTINUE 


GLSRAN2 


COMMENT — NOW HAVE THE PEOUIREO 

S AND R ARRAYS. 

GLSRAN2 

220 

C 




GLSRAN2 




NWNFUNC*(NFUNC+l)/2 


GLSRAN2 




CALL GLSCOVl ( F f S , R # W , B ^ Y^NFUNC » N VV#I CODE# SUM Y# YSQSU M, I ERR ) 

GLSRAN2 




IEND-MMORD+1 


DATM0D2 




IEND-IEND+(IEN0+l)/2 


DATMDD2 

225 



DO 71 I«1»IEND 


DATM0D2 




WRITE (3,789) S(I) 


DATM0D2 



fi' 

CONTINUE 


DATMDD2 


C 




GLSRAN2 




PRINT 201 


GLSRAN2 

230 



NDEG-MORD 


GLSRAN2 




NZM-NOEG+1 


GLSRAN2 




NSM»NZM+2*MnRD 


GLSRAN2 




DO 80 I«1,NFUNC 


GLSRAN2 




IF (INDEX(I).LE.NZM) GO TO 

74 

GLSRAN2 

235 



IF (INDEXd ) .LE.NSM) GO TO 

75 

GLSRAN2 




NOTF-INDEX (I )-NSM 


GLSRAN2 




NOTFl-NOTF+l 


GLSRAN2 




NTG-0 


GLSRAN2 




DO 72 J-1,M0RD 


GLSRAN2 

240 



LDEG«(N0TFl-NTG)/2+J 


GLSRAN2 




IF (LDEG.LE.NDEG) GO TO 73 


GLSRAN2 




NTG-NTG+ (MORD-J i*2 


GLSRAN2 



72 

CONTINUE 


6LSRAN2 



73 

CONTINUE 


6LSRAN2 

245 



LOPD-J 


GLSRAN2 




LE0»M0n(N0TFl,2 ) 


GLSRAN2 




GO TO 76 


GLSRAN2 



74 

CONTINUE 


GLSRAN2 




L0RD«0 


GLSRAN2 

250 



LEO-0 


GLSRAN2 




LDEG«INDEX(I )-l 


GLSRAN2 




GO TO 76 


GLSRAN2 



PKnCRAK RLSSAN2 


7^/74 


nPT«l 


FTN 4.7+485 


80/01/29. 14.11.35 


255 


260 


265 


270 


o 

o 

275 


280 


285 


290 


75 CONTINUE 
NDSF»INDEX(I )-NZM 
LDEG- (NDSF + 1 )/2 
LE0-M0n(N0SF+l,2) 

LORD-LDEG 

76 CONTINUE 

K«I + (I*(I-1) )/2 

WRITE (6»200) e{I)>SC<)»S0RT(S(K))>I,R(I),INDEX(I),LDRD#LDEG/LE0 
80 CONTINUE 
C 

GD TO 95 
90 CONTINUE 

DRINT 110, KC0UNT,JCDUNT,NFUNC,INDEX(NFUNC),I1,I, JJ, J,IVV#IV 
ST0P6 
C 

95 CONTINUE 

DO 900 I«1,NFUNC 
K«=I + (I*(I-1) )/? 

WRITE(3,789) B(1),S(K) 

C WRITE (3,789) B(I),S(K) 

900 CONTINUE 

NNDEG-NFUNC-1-MM0RD*( NMDRD+1) 

VARDATA - (YSOSUf^-SU^Y + SUMY/NDATA) / (NDATA-1) 

BR1»R{ 1 )*B( 1) 

BR«0. 

DO 515 I«1,NFUNC 
BR«PR + 3( I)*R(I) 

515 CONTINUE 
C 
C 

A1-(B(1)*SUMY)/ (NDATA-1 ) 

A?«-(SUMY*SUMY) /(NDATA+{ NDaTA-1) ) 

PRINT 10 05* NPATA,NFUNC »NNDEG,SUNY,SUMY*SUMY, YSQSUM, BR, BR1,A1,A2 

1005 format (♦!♦, ♦fundamental statistical parameters*// 

11X,*T0TAL MEASUREMENTS(NDATA)*,T40,*« *,16/ 

21X,*NUMBER OR MODEL C OEF F I Cl ENTS ( NFUNC ) *, T40, ♦ *= *,14/ 

31X,*DEGREE dp MODEL (NNDEG)*,T40,*« *,13/ 

4iy,*SUMY*,T40, *■ *,E15.8/ 

51X,*SUMY SOUARED(SUMY X SUMY )*, T40, *- *,E15.8/ 
61X**YS0SUM*,T40,*» *,515.8/ 

71X,*YEXPS0SUM(8R)*,T40,*« *,E15.8/ 

81X,*P(1) X R(l) (BR1)**T40,*« *,E15.8/ 


GLSRAN2 
GLSRAN2 
GLSRAN2 
GLSRAN2 
GL3RAN2 
GLSRAN2 
GLSRAN2 
GLSRAN2 
GLSRAN2 
GL3RAN2 
GLSRAN2 
GLSRAN2 
GLSRAN2 
GLSRAN2 
GLSRAN2 
GLSRAN2 
GLSRAN2 
DATMD02 
DATMDD2 
DATMDD2 
GLSRAN2 
DATM0D2 
DATM0D2 
DATM0D2 
DATMOD2 
DATMDD2 
DATMQD2 
DATMDD2 
DATM0D2 
DATMDD2 
DAT MOD 2 
DATMOD2 
0ATM0D2 
0ATMDD2 
DATM0D2 
DATMDD2 
DATMDD2 
DATM0D2 
DATM0D2 
DATM0D2 
DATM0D2 
0ATMDD2 




PROGRAM GLSRAN2 


74/74 OPT-1 


FTN 4.7+405 


80/01/29. 14.11.35 


295 


300 


305 


310 


315 


320 


325 


330 


335 


91X,+Al*>T40#+« ♦,E15.8/ 

♦ lX>+A2=t'#T40»-*'« ♦>E15.8) 

C 

c 

VA'RERR-( YSQSUM-BR) / (NDATA-1) 

VARMOD-VARDATA-VARERR 

C 

HDF-NFUNC-1 

XMSH-VARMQD/MDF 

lERDF-NOATA-NFUNC 

X*1SE«VARERP/IER0F 

C 

C 

COMMENT — CALCULATE THE DEGREE VARIANCES (AVERAGE SQUARE) OF THE 
C SPHERICAL HARMONIC MODEL. 

COMMENT — LET THE FIRST FIVE(5) ELEMENTS OF THE ARRAY R CONTAIN 
C VARDEGCl THROUGH 5)# FOR 5-TH DEGREE MODEL. 

DO 525 I«1#MM0RD 
F(I)»0, 

R(I)»0. 

COMMENT — INN IS THE DEGREE OF THE COEFFICIENT. 

INN«I 

IJ-I+1 

DO 524 J«1>IJ 

COMMENT — JMM IS THE ORDER OF THE COEFFICIENT. 

JMM-J-1 

COMMENT — IF JMM-0, COEFFICIENT IS ZONAL. 

C IF JMM»INN» COEFFICIENT IS SECTORAL. 

C OTHERWISE THE COEFFICIENT IS TESSERAL. 

IF (JMM.EQ.O) GO TO 518 
IF (JMM. EQ. INN) GO TO 520 
C 

COMMENT — CALCULATE B INDEX FOR TESSERAL COEFFICIENTS. 

NT-INN-JMM 

COMMENT — JMMM IS THE NUMBER OP PRECEDING ROWS CONTAINING TESSERAL 
C FUNCTIONS. 

JMMM- JMM-1 

IF (JMMM.EQ.O) GO TO 517 
DO 516 II-I^JMMM 
NT-NT+(MMORD-II ) 

516 CONTINUE 

517 CONTINUE 


DATM0D2 

DATM0D2 

DATM0D2 

daihqd^ 

bATMbb2 
DATM0D2 
DATM0D2 
DATM002 
D ATM 00 2 
DATM0D2 
D ATM 002 
DATM0D2 
DATM0D2 
DATM0D2 
DATM0D2 
DATM0D2 
DAT MOD 2 
DATM0D2 
DATM0D2 
DATM0D2 
DATM0D2 
DATMO^ 
DATM0D2 
DATM0D2 
DAtM0D2' 
DATMDD2 
DATM0D2 
DATM0D2 
DATMbbZ 
DATM0D2 
DATM0D2 
DATM0D2 
DATM0D2 
DATM0D2 
DATM0D2 
DATM0D2 
DATM0D2 
DATM0D2 
DATM0D2 
DATMDD2 
DATM0D2 
DATM0D2 



PROGRAM GLSRAN? 


74/7A OPT«l 


FTN ^.7+^85 


80/01/29, lif.11.35 


o 

ro 


,3A0 


3A5 


350 


355 


360 


365 


370 


375 


KEVEN«3*MMQRD+2*NT 
KQDD-l+KEVEN 
GO TO 521 

5_18_CONTINUE 

C 

COMMENT > CALCULATE B INDEX FOR ZONAL COEFFICIENTS. 

KEVEN-INN+1 
KODD — 99 
GO TO 521 
.. 520, CONTINUE 
C 

COMMENT — CALCULATE 8 INDEX FOR SECTORAL COEFFICIENTS. 
KEVEN-MMDRD+2+JMM 
KODD-KEVEN+1 

521 CONTINUE 
C 

COMMENT — CALCULATE THE "SQUARE ROOT" OF THE EVEN AND ODD TERM 
C CONTRIBUTIONS OF THE DEGREE VARIANCES. 

EVEN«B(KEVEN) 

REVEN-R(KEVEN) 

IF (KODD. EQ. -99) GO TO 522 
ODD-B (KODD) 

RODD-R(KODD) 

GO TO 523 

522 CONTINUE 
□DD«0. 

RODD-0. 

523 CONTINUE 

F(I)-F(I ) + EVEN*REVEN-»-OnD*RODD 
R (I )»R (I )+EVEN* EVEN +000 +000 

PRINT 1001, I, J,KODD»ODO,RODD,KEVEN,EVEN,REVEN,F (I ) 

1001 FORMAT ( 1X,*I-*,I3,+ 

1* KODD-*, 13,* 0D0-*,E15.e, 

2* R0DD-*,E15.e,* KEVEN«*,I3, 

3* EVEN-*,E15.e,* REVEN-*,E15.9, 

A* F(I)-*,E15.6) 

52A CONTINUE 

BR1-8R1 + F ( I ) 

COMMENT — F(I) CONTAINS VALUES FOR THE MEAN SQUARE DUE TO 
C COEFFICIENTS PER DEGREE. 

COMMENT — F(I+NNDEG) CONTAINS VALUES OF DEGREES OF FREEDOM 
C FOR THESE MEAN SQUARE CALCULATIONS. 


DATM0D2 

DATM0D2 

DATMDD2 

DATM0D2 

DATM0D2 

DATM0D2 

DATM0D2 

DATMDD2 

DATM0D2 

DATM0D2 

DATM0D2 

DATM0D2 

DATM0D2 

DATMDD2 

DATMDD2 

0ATM0D2 

DATM0D2 

DATM0D2 

DATM0D2 

DATM0D2 

DATMDD2 

DATMDD2 

0ATM0D2 

DATM0D2 

DATM0D2 

DATM0D2 

DATM0D2 

DATM0D2 

DATM0D2 

DATM002 

0ATM0D2 

DATM0D2 

DATMDD2 

0ATMDD2 

0ATM0D2 

DATM0D2 

DATM002 

DATM0D2 

DATM0D2 

DATM0D2 

DATM0D2 

DATM0D2 



PROGRAM GLSRAN2 




OPT-1 


FTN A.7+A85 


80/01/29. 14.11.35 


380 


385 


390 


395 


POWER PER DEGREE 


o- 

oo 


■400 


405 


410 


415 


420 


F(I)-F(I )/(NDATA-1) 

F(I+NNDEG)»2.+INN+1. 

F{I)-F(I)/F{I+NNDEG) 

COMMENT — S(I) CONTAINS VALUES OF THE ZONAL 
S(I)« B(I+i)*B(I+l)/R(I)*100. 

R(I)»R(I)/(2*INN+1) 

525 CONTINUE 

PRINT 1002> BR#BR1 

1002 FORMAT ( IX, ♦8R F15 . 8, ♦BR !■*» E 1 5 . 8 ) 

TPOWER-0. 

DO 550 I-1,NNDEG 
TPDWER-TPOWER+R (I) 

550 CONTINUE 
C 

PRINT 590, VARDATA,VARMOD, XMSM,MDF,VARERR,XMSE,IERDF, 
1VARM0D/VARDATA,TPCWER,A1,A1+A2, A1+A2 
590 FORMAT ( // /T37# ♦PERC ENT . ♦# T4 9, ♦ VARMOD*, T66»*DE G. CONTRIB.*, 
*T83,*ACCU.MULATIVE*,T100,*MEAN SQUARE*, T 117, *DE G . FREEDOM*/ 
♦T37,*Z0N. POWER*, T49,*DEG. CONTRI B .*, T69, ♦+ A2*/ 

♦1X,*VARDATA« *,E15.8/ 

11X,*VARMDD- *,E15.8,T98,E15 .8,T117, 15/ 

21X,*VARERR» +, E15.8,T98, E15. 8,T117, 15/ 

31X,*RATI0 - *,E15.8/ 

41X,*TP0WER- *,E15.8// 

5T43,*A1«*, E15.8,T64, E15.8,T81, E15.8) 

C 

ACCUM«A1+A2 

DO 535 I-1,NNDEG 

CAPPA-Fa)*F(I+NNDEG) 

ArCUM-ACCUM+CAPPA 

PRINT 595, I,R(I),S(I),CAPPA,CAPPA+A2,ACCUM,F(I),F(I+NNDEG) 

595 FORMAT ( 1 X, * VARDE G ( ♦, 1 2, * ) - *, E 1 5 . 8, T36, F 10, 5, T47, El 5 . 8, T64, E15 . 
1TB1,E15.8,T98,E15.8,T117,F5.0) 

COMMENT — LET R(I) NOW CONTAIN PERCENTAGE POWER, 


535 


585 

545 

950 

999 


P(I)*(R(I)/TPDWER)*100. 

CONTINUE 

DO 545 I*1,NNDEG 

PRINT 585, I,R(I) 

F0RMAT(1X,I2,5X,*PERCENT 

CONTINUE 

CONTINUE 

CONTINUE 


POWER- *,F10.5) 


DATM0D2 
DATM0D2 
DATM0D2 
DATM0D2 
DATMDD2 
DATM0D2 
DATM0D2 
DATH0D2 
DATMDD2 
DATM0D2 
DATM0D2 
DATMDD2 
DATM0D2 
DATM0D2 
DATM0D2 
DATM_0D2 
DATMOD’2 
DATM0D2 
DATM0D2 
0ATM0b2 
DATM0D2 
DATM0D2 
DATMOD‘2 
DATM0D2 
DATM0D2 
DATM0D2 
OAT MOD 2 
DATM0D_2 
DATM0D2 
DATMDD2 
DATM0D2 
6,DATM0D2 
DATM0D2 
DATM0D2 
DATMDD2 
DATM0D2 
DATM0D2 
DATM0D2 
DATM0D2 
DATM0D2 
GLSRAN2 
GLSRAN2 



PR3G»?A^! GLSRAN2 


7A/74 


DPT-1 


FTN A.7+A85 


80/01/29. 19.11.35 


925 


930 


935 


990 


100 FORMAT 
1E15.8) 
105 FORMAT 
110 FORMAT 


STOP GLSRAN2 

C GLSRAN2 

C GLSRAN2 

(lX/tZl» E15.8#5X#+22- E15 . 8» 5 X# ♦ Y- ♦> E15 . 8> 5X» *ER» ** GLSRAN2 

GLSRAN2 

(T10#+T- +/E15.6) GLSRAN2 

(♦1*#*ST0P6 INDICATES A POTENTIAL RUN-AWAY LOOP SITUATION EGLSRAN2 
IXISTS IN LOOP - 70 IN GLSR AN2 .♦ /I X, *THE FOLLOWING PARAMETERS ARE PGLSRAN2 
PRINTED AS DIAGNOSTIC AIDS — KCOUNT# J COUNT » NFUNC^ INDE X ( NPUNC ) ^ I 1 / IGLSRAN2 
3/ JJ# J#IVV>IV*/1X,10(I5^5X) ) GLSRAN2 

175 format (♦l + »//////////T25,*BEGIN PRINT FOP CASE NUMBER 1 3 / / / / / / / G LSRAN2 
1///) GLSRAN2 

200 FORMAT (TIO, El 5 . 8# T30# E 1 5 . 8# T50> E 15 . 8, T2, 1 3, T72> E 1 5 . 8, T90> I 9, TlOl, GL S RAN 2 

1I2»T110,I2,T120,I1) GLSRAN2 

201 FORMAT (//T19»*EST. C OE F . T39 » *8 ET A VAR . ♦ » T 53 # T . DEVIATION*, GLSRAN2 

lT75,*P-VECT0R*,T90,*INDFV*,T99,*0RDER»,T10e,*DEGREE*,T118, GLSRAN2 

2*EVEN-0*/ GLSRAN2 

3T17,*B *,T39,*C0VAR (1 ,I ) *, T53, *50 RT (C 0 V AR ( I, I ) )*, T90» *APRA Y*, T119, GLSRAN2 
9*rDD«l*/) GLSRAN2 


215 FORMAT (///1X>A7) 
789 F0RMAT(2 (5X,E15.8) ) 
END 


GLSRAN2 

DATM0D2 

GLSRAN2 


o 

-p> 



SYMBOLIC 

REFERENCE 

MAP (R-1) 

TPY 

POINTS 



502 

GLSRAN2 



5TARLFS SN 

TYRE 

RELOCATION 

356 

ACCOM 

REAL 


333 

A2 

real 


661 

BETA 

REAL 

ARRAY 

3 30 

BRl 

REAL 


A 

OX 

REAL 

/ / 

3‘^1 

EVEN 

real 


276 

I 

INTEGER 


263 

ICQDE 

INTEGER 


252 

IDATA 

INTEGER 



2C332 

A1 

REAL 


55137 

B 

REAL 

ARRAY 

2G331 

BE 

REAL 


20357 

CAPPA 

REAL 


6 

ER 

PEAL 


2C360 

F 

real 

ARRAY 

2C256 

iCASE 

INTEGER 


20275 

ICOUNT 

integer 


20313 

LEND 

INTEGER 




SOBROUTiNE GLSCORl 


7A/7-<* OPT-1 


PIN 4. 7+^85 


80/01/23. 19.03.20 


1 


5 


10 


15 


20 


o 

cn 


25 


30 


35 


AO 


SUBROUTINE GLSCORl { F^S>R>Wi Y# NV# 1C ODE# S UMY# YSQSUM. I ERR ) 
DlhFNSICN P(N),S(NV)>R(N)^B(N) 

C 

COMMENT -- INITIALIZATION LOOP FOR THE R^B> AND S ARRAYS. 

CO 25 I-l^N 
R(I )«0. 

B(I)»0. 

25 CONTINUE 

DC 35 I-1,NV 
S(I )-0. 

35 CONTINUE 
SUE* Y-0. 

YSQSUM-0. 

RETURN 

C 

c******* ♦♦♦♦ 

c 

ENTRY GLSSLMl 

COHN ENT — SUMMATION OF THE R AND S ARRAYS. 

IF (ICODt.ECi.O) GO TO 150 
DC 125 I-liN 
R(1)»R(I)+F(1)*W^Y 
125 CONTINUE 
C 
C 

COMMENT — CALCULATE THE SUM OF IHt Y»'S AND Y SQUARED SUMMED. 
SUNY-SUKY+Y 
YSOSUM«YSQSUM+Y*Y 
150 CONTINUE 

COMMENT — CALCULaTl THE UPPER FULl TRIANGLF OF THE S-MATRIX AND STORE 
C THE kESULT in THE ON E-Di MENS 1 ONaL ARRAY> S> DIMENSIONED 

C BY NV — SEE PARAMETER LliT IN GlSRANZ. 

DO 175 J«1^N 
DO 175 1«1,0 
K-J + U-M J-1) ) /2 
S(K)*S(K)+F(I)*P(JJ*W 
175 continue 
RETURN 
C 

C* + ***** *♦♦♦ 

c******* ♦♦♦♦ 



SUcKL-UlNL GLSCQki 


7^/7A 0FT*1 


FTN 4,7+^d5 


3C/01/23, 19*03*20 




50 


55 


60 


o 

CTl 


65 


C 

LNTFY GLSCLVi 

CCMkcM — find INVLRSt uF S AND STukt IN S. 

C FIND ESTlhATICN PAkAhfcT£k6> c»SP WHhRE S IS NOW THE 

C CCVARIANCF MATRIX OF PARAMETRIC ESTIMATION. 

iCP-I 

CALI SFCIMf (N,S,JCP>DhUISCALE#IEkR) 

IF (ICDUE.EO.O) GO TO 2f0 
DL 225 I«1^K 
E (I )»0. 

CL 225 J*i^N 

1L*1 

JD- J 

IF (i .LF.J ) GO TO 200 
1 Ci^J 
Jij» i 

2l 0 CGNTINLt 

K*iO+( J0*( Ju-1 ) )/Z 
t(i )«B(I)+i(K)4N( J) 

225 CL^Tj.NUF 
250 COI.TISUF 
RiTUFN 
INC 



S YMbDLiC 

RtFl PtNtL 

MAP (R-1) 


NTkY 

P 0 } f T S 




3 

GlSCCRl 

101 

GLSCDVi 


AKiABLES SN 

TYPE 

RtUKATION 

0 

B 

U £ A 1. 

ARF AY 

r.p. 

0 

r 

RtA L 

ARRAY 

F.P. 

t 

i c c u t 

IN 1 t; G i ► 


f 

167 

i r 

iNTt GEK 



166 

ISCmLF 

IM EC-5R 



170 

J[ 

INTEGER 



0 

i\ 

I NT EGER 


F • F • 

U 

p 

HP A L 

ARRAY 

F.P. 


GLSSUMl 


io5 

uET 

RfaL 



161 

I 

INTEGER 



0 

I ERR 

INTEGER 


R.R. 

16A 

i, uP 

INTEGER 



162 

J 

INTEGER 



163 

K 

INTEGER 



C 

NV 

INTEGER 


F.P. 

0 


REAL 

ARRAY 

F.P. 



iUBROUTiNE SIhDATl 


74/7A uPT-1 


FTN A.7+A85 


60/01/23. 19.03.20 


1 


5 


10 


15 


20 


o 

--j- 


25 


30 


35 


Hit 


SUBROUTINE SIK.DATl (BETAjNiF) 

COMKCN PI,Xl>X2,YVAk,DX#X>ER> IFUNC>Y,Z1,Z2> Z>M0RD 
DIFENSIDN BITA(N),F(N) 

CALL SIRDAT2 

CQMFtNT — SEE REALDAT FOP EXPLANATION OF DATA TO BE SIMULATtD AS INPUT. 
IF (MORD.NE.O) CALL SIMDAT3 

IF (IFUNC.EQ.l) GO TO 25 

IF (1FUNC.EQ.2) GO TO 125 

IF (I FUNC.EC.A ) GO TO 250 

IFL'NC-909 
RETURN 
C 
C 

CDMFENT — CALCULATE THE LINfcAK POLYNOMIAL FUNCTIONS OF DEGREE •’NP«, 

25 CONTINUE 

CALL POLY (F>bETA>N) 

50 CDMTINUF 
Y»Y+FR 
RETURN 
C 
C 

COMMENT — CALCULATE THE FIRST •'NP« ZONAL SPHERICAL HARMONIC FUNCTIDNS. 
125 CONTINUE 

CALL SPHARN (F>N) 

135 CONTINUE 

Y-0. 

DO 150 1«1>N 
Y*Y+BETA(1J*F{1) 

150 CONTINUE 
GO TC 50 

COMMENT — END SPHERICAL HARMONIC CALCULATIONS. 

C 


C 

COMMENT — CALCULATE THE "NP" FOURltR rONCTI0NS» OTHER THAN F(l)»l, 
250 CONTiNUr 

CALL FDRFNCl (F»N) 

GC TO 135 

COMMENT — END FDURIEk CALCULATIONS. 

C 

C 


END 




U6KDUUNE RlALb/^T 

L'?T«1 


FTN 

A.7+A85 3C/D1/23. 19.03. 20 

1 


sue, cnuTiKc 

FFALDaT (F,N) 






CCrhLN Pl>Xl,X2»YVAR,DX>X#tK, 

IFUNC* Y,21»Z2*2*N0RD 



L 

C 

LIPEKiILN f(N) 




t> 







Cl 

jMfllM — EXPlaKATIOh JF INPUT. 





L 


X 

Z 


Y 


C 

lFUNC«i> 

INDtFF.NDFNT 

NON? 


DEPENDENT 


c 

1 


VARi ABL E. 



VARIABLE 

X tJ 

w 

r 

1 Fire «2, 

CU-LAl 1 TUCt 

LONGITUDE 


DLPENDFia 


c 


(RADIANS ) 

(RADIANS) 


VARIABLE 


c 

(• 





(OZONE) 

15 

c 

i FLf C-3> 

LATITUD'^ 

NONF 


DEPENDENT 


c 


(RaDImNS) 



DEPENDENT 


c 





VARIABLE 


c 





(3UV-GRIDED 


c 





MODEL DATA 

20 

c 





CORRECTED 

108 

c 

c 





FOR DOBSON) 


c 

i FUNC*4# 

SCALfcU FOURIER 

mONE 


APPROPRIATE 


c 


AnGLE(RACIANS » 



DEPENDENT 

25 

c 

c 





VARIABLE 


c 

iF (MGRD.EC 

.0) GO TO 15 






Rt*r. (2) X,Z>Y 




3 j 


IF ( f a F ( ? ) ) 

2 5 » 5 n 






15 CLi'TINL' 







pr/r!(2) x,Y 






r 

If (faM2)) 

25*50 




35 

c 

2 5 CcMlNUl- 







lFUNt«'>»9e 






f 

nF TUFN 





^0 

V 

c 

5 - Cl I'T J wur 







IF ( Y F Cr* C « F 
1 * 

O.i) GT To 75 






SUBKOUTlNi: REALDAT 


74/7^ OPT-1 


FTN ^.7+485 


eO/01/23. 19.03^20 


45 


50 


55 


60 



70 


75 


IF (1FUNC.E0.2) GO TO 150 

IF (1FUNC.EQ.3) GO TO 250 

IF (JFUnC.EO.4) go to 350 

IFUNC-995 
RETURN 
C 
C 

CCMhENT — CALCULATE LINEAR POLYNUrtlAL FUNCTIONS THROUGH DEGREE »NP«-. 
73 CCNTINUE 

DU 125 I-2#N 
F(I )«F(1~1)*X 
125 CONTINUE 
C 

RETURN 

C 

C 

COMMENT — CALCULATE THE "NP” SPHFkICaL HARMONIC FUNCTIONS. 

150 CONTINUE 

CALL SPHARM (F, H) 

RETURN 

C 

C 

COMMENT — CALCULATE F ( 2 ) -COS ( 2 + LAT ) , WHERE X-LAT. 

250 CONTINUE 

F (2)«CCS (2*X ) 

RETURN 

C 

A 

V 

COMMENT — CALCULATE THE »NP" FOURIER FUNCTIONS^ OTHER THAN F(l)-1. 

350 cor TiNue 

CALL FDRFNCl (F,Nl 
RETURN 
C 
C 

END 


SYMBuLiC REFERENCE MAP (R-l) 



SUaKDUTlht POLY 


74/7A CPT»1 


FTN 4.7+A85 


80/01/23. 19.03.20 


1 SbBFDUTIhE POLY (F^A>N) 

CDP.MCN Pl»Xl#X2>YVAF,DX#X»Ek#lFUNC>G 
DIMENSION F(N)/A(N) 

NP-N-1 

5 G-A(NP+1) 

If (NP.EQ.O) RETURN 
DC. 15 M-1,NP 
J-KP+l-M 
F(M+1)«F(M)*X 

- J.0 G-A(J)+X»G 

15 CONTINUE 
RETURN 
END 



SYMBOLIC 

REFERENCE 

MAP (R«l) 

TRY 

POINTS 



-.3. 

POLY 



RlAbLES SN 

TYPE 

RELOCATION 

0 

A 

real 

ARRAY F.P. 

6 

Ek 

REAL 

/ / 

10 

6 

REAL 

/ / 

31 

J 

INTEGER 


0 

N 

INTEGER 

F.P. 

0 

PI 

REAL 

/ / 

1 

XI 

REAL 

/ / 

3 

YVAR 

REAL 

/ / 

ATEM 

SNT LABELS 



U 

15 



)UPS 

LABEL 

INDEX 

FROM-TO LENGTH 

20 

15 

M 

7 11 5B 

IMMUN 

BLOCKS 

LENGTH 



/ / 

9 



A 

DX 

REAL 


/ / 

0 

F 

real 

ARRAY 

F. 

7 

IFUNC 

INTEGER 


/ / 

30 

M 

INTEGER 



27 

NP 

INTEGER 



5 

X 

REAL 


/ / 

2 

X2 

REAL 


/ / 


PROPERTIES 

INSTACK 



SUBROUTINE SPHARN 


7A/7A OPT-1 


FIN A.7+AB5 


80/01/2S. 19*03*20 


1 


5 


10 


15 


20 


25 


30 


35 


^0 


SUBROUTINE SPHARtt (F/N) 

CDNMDN PIjX1^X2> YVAK#DX^ X^ER>1FUNCjY,Z 1*Z2>Z^M0RD 
DIRENSION F(l) 

DINENSIDN Q(13>13) 

C 

C 

C ♦ ^ ^ 4> 

comment — CALCULATE THE SECOND THROUGH THE (NDEG+1) - TH F-FUNCTIDNS. 
C 

MPLUS-MDRD+1 

r.F»N-i 

NDLG«NP“MDRL*MPLUS 

CALL LEGNDR4 ( M PLUS > NDEG , X , 0 ) 

M»0 

00 25 1-l^NDEG 
FlI+i)»Qll#l+l) 

25 CONTINUE 

C*** + + *** 

C 

c 

c >» *♦♦♦»♦♦ 

COMMENT — CALCULATE THE 2*MQP.D SECTORAL FUNCTIONS. 

C STORE results IN THE F ARRAY# ELEMENTS (NDEG+2) THROUGH 

C (NDEG+1+2+M0RD) . 

C 

NM-NDEG+3 

NN2»2*MURD+NNl-2 

M-C 

DO 50 I«NN1#NN2#2 
M + l 

FS-SCkT(2./FFAC(M + M) ) 

F(l-l)«FN*CuS (M*Z) 

F(I )«FN*SIN(M*Z) 

50 CONTINUE 

IF (MORD.EO.l) RETURN 

C***»***» 

C 

c 

C*****^t* 

comment — CALCULATE THE NUMTtS-N-NN2 TESSERAL FUNCTIONS. 

C 



SLI6RCLTINL SFhARM 


7^/74 CFT-1 


NM-NN2 + 2 
NK2-N 

45 M»C 

NN-NDtG 

Dc 75 l*KNl>NK2>2 
IF (NN.LT.NDEG) GO TO 70 

M-F.41 

50 NN«F 

70 CDMlNUt 
Nh»NN4l 

fS»SORT(2.*RFAG(NN-N) /RrAC(NN+M) ) 
FN*Q (M+1#NN+1) *FS 

55 F ( l-l)*=FN*CiJS {rt*Z) 

FU )«FN»SiN(M4Z) 

75 CONTINUE 
KETURN 
FKC 


ro 


SYMBOLIC RtFFFENCE MAP (R«l) 

nTRy points 

3 SFHaRR 


ARIA6LE5 

SN TYPE 

RLLDCATICK 


4 

DX 

REAL 

/ / 

b 

0 

F 

REAL 

ARRAY F,P, 

153 

152 

F5 

REAL 


147 

7 

J func 

] M EGEP 

/ / 

inb 

14 

ilURu 

INT! GtN 

/ / 

143 

0 

H 

I N T 1 6 E R 

? . P • 

145 

154 

NN 

INTEGER 


150 

151 

NN2 

INTi.GEF 


144 

0 

PI 

real 

/ / 

155 

5 

X 

REaL 

/ / 

J, 

2 

X2 

<<EAL 

/ / 

10 

3 

yvak 

real 

/ / 

13 

11 

11 

PtAl 

/ / 

12 


FTN 4,7+465 


60/01/23. 19.03.20 


EP 

REAL 


/ / 

FN 

REAL 



I 

INTEGER 



M 

INTEGER 



MPl US 

INTEGER 



NO E G 

INTEGER 



Uhl 

INTEGER 



HP 

Integer 



Q 

real 

ARRAY 


XI 

REAL 


/ / 

Y 

REAL 


/ / 

Z 

REAL 


/ / 

Z2 

REAL 


/ / 



fUNCTIDN Rl-AC 


74/7^ UPT-1 


FTN 4.7*48i> 


80/01/23. 19.03.2C 


1 


5 


LLiLLLc function kFAC (Hu) 
PFaC«1, 

1KN0.lt. 2) GO TO 11 
UD 9 I-1,N0 
KFAC«RFAC*I 
9 CLMJNUL 
11 FtTUkN 
fcND 


SYMBOLIC PtFEFtNCE MAP (R-1) 

tky points 

5 RfAC 

RiABLtS SN TYfE RtLOCATION 

30 1 INTEGER U NO INTEGER F.P. 

26 RFAC UDUfiU 

ATfchfc'NT LAuELS 


0 9 


24 

11 

OPS L/BEL INCEX 

FROM-TO 

LENGTH 

PROPERTIES 

16 9 1 

4 6 

58 

INSTACK 

ATISTICS 
PROGRAM length 

31B 

25 



620008 CM USED 



SUBROUTINE fORFNCl 


7 ^/ 7 <. 


CJPT-1 


FTN 4,7+4o5 


8C/01/£3. 19#03.20 


1 


0 


10 


15 


20 


25 


SUBROUTINE FOftENCl (F#N) 

COMMON n,Xl>X2>YVAR,DX,X 
DIMENSION F(N) 

C ******** ****i‘>r******i‘***** ********************* ******* 


C 

COMMENT — NP»N-1 FOURIER FUNCTIONS ARE CALCULATED PER ♦ 
C CALL TO subroutine FDRFNCl* ♦ 
C — THESE NP FUNCTIONS ARE OF THE FORM, ♦ 
C F(2*1)«C0S (I*X) ♦ 
C AND ♦ 
C F(2^1+l)» SIN ♦ 
C FOR 1-1, M ♦ 
C i<HtRE rt»NP/2* * 
C 


Q *4.^ ****** *****^***********m****** ****** *************** 

C 

C 

M-(N-l)/2 
DO 25 I«1,M 
F(2*1)“CDS(I*X) 

F (Z + l+D-SINd + X) 

25 CDKTlNUt 
C 

RLTURN 

END 


SYMBOLIC RtFEFENCF MAP (R»l) 

,TRY POINTS 
3 FOkFNCl 


^RIABLES SN TYFE RELOCATION 


A 

DX 

REAL 

/ ! 

C 

F 

F£AL 

ARRAY 

F.P, 

2A 

I 

INTlGEF 


23 

M 

INTEGER 



0 

N 

integer 

F.P. 

0 

PI 

REAL 


/ / 

5 

X 

REmL 

/ / 

1 

XI 

real 


/ / 

2 

X2 

REAL 

/ / 

3 

YVAR 

REAL 


/ / 



SUBRubllNE 

1 

5 

10 

15 

20 

25 

30 

35 

^0 


LnGNOkA 74/7-% C?T»1 FTN 4*7+485 80/-01/23# 


SUBROUTINE LEGNDR4 ( NORD# NOE G#CQLAT> P ) 

DIKENSIDN PI13j1) 

DDUELE Q113#13) 

. DOUBLE XjilKE 
MQkD«NORD-l 
Q(l> 1)-1* 

F(l>l)-Q(l^l) 

IF (NDEG.EQ.O) RETURN 
X«CDS (CDLAT) 

SlNE-DSQkl (1*-X+X) 

0(1,2 )«X 
P(l,2)-Q(l>2) 

IF (NDLG.tQ.l*AND,MDRD.cQ«0) RETURN 
N«C 

50 CONTINUE 

If lflOKD.Nt.OJ GO TD 150 
N«N+1 

75 CONTINUE 

CONSENT — CALCULATt ZERO ORDER TERN OF DEGREE N+1 WITH THE TWO 
C PREVIOUS ZERO ORDER TEPfiS. 

Q(l,N+2)«( {2+N+l)+X+Q(l,N+l)-N*Qll,N> )/(N+l) 
P(l,N+2)-0ll,N+2) 

IF (HORD,EO,O.AND.NOEG.£Q.N+1 ) RETURN 
GO TO 50 
150 CONTINUE 
N"N + 1 
M«C 

225 CONTINUE 

COMNENT — CALCULATE HIGHt.K THAN ZERO ORDER TERMS OF DEGREE N. 
M-N + 1 

IF (X,tQ.l.,aR.X,EO.-i. ) GO TC 2 p0 
TO'O. 

IF {Im*GF,M + 2) TQ«Q(M + 1,N-1) 

Q(N+l,N+l) «T0+ (2+N-i)+£lNE*U ( M,N) 

GO ID 30( 

2 50 CONTINUE 

0(N+i,N+l)»C. 

300 CONTINUE 

P(N+1,N+1)-Q(M+1,N+1) 

IF (M.EQ.hLRO.ANO.N.tQ.NDSG) RETURN 
IF (N.FQ.N.OR.M.EQ.MORD) GO TO 75 
GD TD 225 


19 * 43^20 



iUBKLUTlNL lEG:^uR^ 


74/7^ CiPT-1 


FTN 4.7+485 


80/01/23. 19.03.2C 


EM 


KD Nk. SEVkuTY DETAILS DlAGhQjlS D+ PROBLEM 

12 I P ARRAY RLFFhENCE OUTSIDE DlflENSlON BOUNDS. 


SYMbDLlC PEFERF.NCL MAP (R-1) 

ilTJtY POINTS 
3 LEGNDK4 


ARIABLES 

SN TYFE 

MLDCATIOi'l 





0 

CULaT 

REAL 

F.P. 

lb3 

ri 

INTEGER 


161 

NOKD 

INTEGER 


162 

N 

INTEGER 


0 

NDcG 

INTEGER 

F.P. 

0 

NDRD 

INTcGER 


0 

P 

REAL 

array F.P. 

165 

0 

C0U9LE 

ARRAY 

157 

SINE 

DOUBLE 


164 

TC 

REAL 


155 

X 

DOUBLE 






XTERNALS 

TYPt 

ARGS 






CDS 

REAL 

1 LIBRARY 


OSQRT 

DOUBLE 

1 LIBRARY 


TATEMENT LABELS 
35 50 

101 225 


40 75 

132 250 


76 150 

136 300 


JATISTICS 

PKOGkAM LENGTH 707B 455 

520006 CM USED 



iUBicDUTINE SUBSl 


7^/7A QPT-1 


FTN 4.7+4B5 


80701723. 19..03.20 


1 


5 


. 10 


15 


20 


25 


30 


35 


4C 


SUBROUTINE SJBSl ( INOEX#Nh>MNORD«N>KORO ) 

DihENSIDN INDEX (NN) 

CDhhENT — DEFINE REQUIRED PARAMETERS. 

NND E G-NN-MMDRD* ( ilMDkD+1 1 -1 

NSZQK-NNDEG+1 

NZS-NSZQN 

ISS«1+N2S 

NSSEC«2»MhDRD 

NS5-NZS+NS5EC 

ITS-l+NSS 

NSTcS-Nh-thSZaN+NSSFCJ 

NTS-NSS+NSTES 

NDE&*>N-MaRC*(MDRD + l)-l 

NF.ZON-NDEG+1 

NZM»NMZDN 

NMSEC-2+MaRD 

NSM«NZM+NMS£C 

COMMENT — CALCULATE THE ELlMcNTS OF THE INDEX ARRAY. 

DL 50 KZDN-1>NZS 
INDEXIKZUN)«KZDN 
50 CONTINUE 

K*KZM . 

DO 100 KSEC-ISS>NSS 
K«K + 1 

iND£X(K5cC)«K 
ICO CONTINUE 

IF (NSTES.EQ.O) GO TO 200 

COMMENT — IRDWS IS THE TOTAL NUMBER OF kOWS OF TESSERAL FUNCTIONS 
C IN THE SUBSET MODEL. 

IRLWS-MMDRD-1 

K«NSM 

CDMMtNT — LEFT IS THl NUMBER DF 'MASTER' SPHERICAL HARMONIC TESSERAL 
C FUNCTIONS REMAINING IN ROW IftOW. 

LEFT-0 

KTES-NSS 

DC 175 1ROW-1>1ROWS 
M-K+LEFI 

COMMENT — NDMTk 15 THE NUMBER OF 'MmSTER' MODEL TESSERALS IN ROW IROW. 

C NOMMTF IS THE NUMBER OF 'SUBSET' MODEL TESSFRALS IN ROW IROW. 

HCMTR-(MaRD“IPnW)*2 
KQhMTR«{MMORD-I ROWJ-2 

DL 150 KSTEP-l^NOMMTR 



SUBKOljTlhE SUBSl 74/7^ DPT-1 


KTES-KTES+1 
K-K + 1 

45 . INDEXCKTES )-K 

L5.0.CCf;TINUE 

LEFT-NDMTK-NDMMTR 
175 CONTINUE 
200 CONTINUE 

50 RETURN 

END 



SYMBOLIC 

REFERENCE 

MAP (R»l) 


URY 

"poi Nl s 




3 

SUBSl 




^RIABLES SN 

TYPE 

RELOCATION 


0 

INDEX 

INTEGtP 

ARRAY F.P. 

125 

-122 

_1R0WS 

. INTEGER 


1J4 

107 

ns 

INTEGER 


120 

121 

KSbC 

INTEGER 


130 

124 

KTES 

INTEGER 


117 

123 

LEFT 

INTEGER 


0 

0 

MORE 

INTEGER 

F.P. 

0 

-112 

. NDEt 

INTEGER 


115 

113 

NMZCN 

INTEGER 


0 

101 

NKUEG 

INT EGER 


127 

126 

NUMTF 

INTEGER 


116 

106 

NSS 

INTEGER 


105 

110 

NSTES 

INTEGER 


102 

111 

NTS 

INTEGER 


114 

103 

N2S 

INI EGER 




lATErtENT LABELS 
0 50 

0 175 


0 100 
ICC 200 


FTN 4*7+4fi5 


BC/01/23. 19.03.20 


IRQW 

INTEGER 


ISS 

INTEGER 


K 

INTEGER 


KSTFP 

INTEGER 


KZON 

INTEGER 


MMQRD 

INTEGER 

F.P 

N 

INTEGER 

F.P 

NMSEC 

INTEGER 


NN 

INTEGER 

F.P 

NDMMTR 

INTEGER 


NSM 

INTEGER 


NSStC 

INTEGER 


NSZUN 

INTcGER 


NZM 

INTEGER 



0 150 



StBkDUTIKE SbBS2 


74/74 DPT«1 


FTN 4,7+4b5 


&0/01Y23, 19*03^20 


1 SUBROUTINE SUBS2 ( INDEX, MQRD) 

COMMENT — GIVEN A NUMBER, M, Of INDEPENDENT OBSERVATIONS, SUBROUTINE 
C SU8S2 CALCULATES THE ARRAT INDEX WHICH CONNECTS THE INDICES 

C Of THE S-MATRICES(UPPER FULL TRIANGLES ONLY) OF THE MASTER 

5 C MODELIORDER-MORO) AND THE SUBSET MDDELISIZE TO BE 

C DETERMINED IN THIS SUBROUTINE BASED ON M)* 

C 

COMMENT — NOTE M MUSI EE .Gt, 1. 

C 

10 DIMENSION INDEX(M) 

NMS EC*2*MCJRL' 

IMAX-MORD+1 

COMMENT — NMSEC IS THE MAXIMUM NUMBER OF SECTORAL FUNCTIONS IN THE 

C MASTER MODEL. 

15 C IMAX IS THE MAXIMUM NUMBER OF ZONAL FUNCTIONS IN THE , . 

. C - MASTER MODEL — REFERRED TO AS NMZON IN SUBROUTINE SUBSul* 

IMA XSO-IMAX+IMAX 

COMMENT — IHAXSQ IS THE MAXIMUM NUMBER OF FUNCTIONS AVAILABLE IN THE 
C MASTER MODEL. 

20 IF (M-lMAXSu) 15,225,250 

15 CONTINUE 

DO 25 1-2, IMAX 

^ If (I+l.GT.M) GO TO 50 

25 CCNTINUE 

25 COMMENT — THIS LOOP SHOULD NUT FINISH NORMALLY. IF IT DOES A 

C DIAGNOSTIC WILL BE PRINTED AND EXECUTION STOPPED. 

PRINT 925 

925 FORMAT (♦!♦,♦ EXECUTI ON STOPPtU IN SUBROUTINE SUBS2*) 

ST0P4 

30 C 

C 

50 CCNTINUE 

NlSSLC«2*(I-2) 

35 ^ISZON-I-l 

COMMENT — NI5SEC AND NISZON ARE THE NUMBER OF SECTORAL AND ZONAL 
C FUNCTIONS, RESPECTIVELY, IN THE INITIAL SUBSET MODEL OF 

C ORDER- 1-2. 

MtIFf -M-NISZON+NISZQN 

4C COMMENT — MUIFF IS THE DIFFERcNCt BETWEEN M AND THE NUMBER OF FUNCTIONS 

C Ih A MCiUtL OF OKDEk-1-2, WHICH CONTAINS (I-1)*(I-1) FUNCTIONS 

MSCC-MDlFF/2 




FTN ^•7+435> 


iUBkOUTlNE SUBS2 7A/7A OPT»l 


8C/01/23# 19*03#20 


^5 


. i>0 




60 


65 

ro 

o 


. 70 


75 


BO 


«ZXiN»l1DlFF-hSEC 
MTES-O . 

CLMEENT — - hS£C AND -NZON ARE THE NUMBER OF EXTRA SECTORAL ANO ZONAL 
C FUKCTIHNS, RESPECTIVELY^ MT£;> IS THE NUMBER OF EXTRA 

C TESSEFAL FUNCTIONS RtUUIREOi IF ANY. 

ITEST-NISSEC+MSEC-NMSEC 
iF UTEST) 75^75^70 
70 CONTINUE 

COMMENT — MSEC+KiSStC IS LARGER THAN NHSEC^ THE MAXIMUM NUMBER OF 
C SECTORAL FUNCTIONS AVAILABLE AN THE MASTER MODEL. 

MStC-NMSEC-NlSSrC 
MTtS-JTEST 
75 CQMINUE 

ITEST ■MZCN+NlSZQi'i-IMAX 
IF (ITESTI 65,65,60 
. to CONTI NLE 

MZDN-IMAX-NISZDiS 
MTtS«MTES+ITEST 
F5 CCNTINUF 

COMrtNT — THE CALCULATION uF MTtS IS NOT ACTUALLY REOUiRED IN ORDER 
C TC FIND NSTES — HOWEVtR, AS A DEBUGGING AID IT IS A 

C USEFUL parameter. 

NZS-KISZDN+MZON 

lSS-NZS+1 

NSS-NZS4MSSEC + MSEC 
ITS-NSS-H 
NS5ct-NS5-NZS 
NSTEl-M-NSSEC-NZS 
NTS -NSS+NSTES 
C 

C * ’r- ♦ ^ * ♦ ♦ * ♦ ♦ ♦ 

c » 

c 

DC 125 KZDN-i,NZS 
INDEX (KZON )»KZON 
1?5 CCNTiNUr 

ir (NSSE( .n.J ) GU TO 155 
K-IMAA 

DC 150 KScC«iSS»NS5 
K-h + 1 

lMtX(KSEC)-K 
150 CuMlNUr 




SUBROUTINE SUBS2 


7^/7A DPT-1 


FTN ^,7+485 


80/01723* 19.T13^20 


85 


90 


95 


100 


105 


liO 


115 


120 


125 


155 CONTINUE 

IT (NST£S.UT*1) GO TO 180 

COhMcNT — iPOWS 13 THE TOTAL NUMBER OF RDkJS OF TESSERAL FUNCTIONS 
C IN THE SUBSlT MODEL. 

IkOWS-I-2 

C NOTE — RECALL THAT NlSZON-I-1 

IF (MSZQN^MSZON.EQ.M) IROWS «IR0WS-1 

K-3+M0RD41 

COMMENT — LEFT IS THE NUMBER OF 'MASTER* SPHERICAL HARMONIC TESSERAL 
C FUNCTIONS REMAINING IN ROW IROW. 

LEFT-0 

KTES-NSS 

DC 175 lROk-1, IROWS 
K-K+LEFT 

COMMENT — NUMTR IS THE NUMBER OF 'MASTER* MODEL TESSERALS IN ROW IROW. 

C NOMMTR IS THE NUMBcR OF 'SUBSET* MODEL TESSERALS IN ROW.-LROW^- 

NOMTP- (M0RL-Ik0W)*2 
NCMMTR-(NiSZDN-IP.DW )*2 

IF (NISZDN*NISZON.£U.M) NuMMTP.-NQMMTR-2 
DO IfcC KSTEP-1>NDMMTR 
KTtS-KTES+1 

IF (KTcS.GT.tt} GO TO 180 — - - 

K-K+1 

1NDEX(KTES)-K 
160 CONTINUE 

LtFT-NDMTR-NaHMTR 
175 CONTINUE 

180 CONTINUE 

RtTUPN 

C 

miZS COMTIMOE 

PRINT 95C'j M 

950 FORMAT (/////1X#*M» IS THE TOTAL NUMBER OF FUNCTIONS AVAIL 

IABlC in THf SPECIFIED MASTER MOOS L. ♦ /IX , ♦THERE P3RE # THE UPPER FULL 
2TRJANGLE OF THt MASTcR MCDcL WILL BE READ DIRECTLY FROM TaPE AND 
3U5tD IN THF FOLLOWING C AL CUL ATI QNS .♦ ) 

GO TO 30U 
250 CCNTINUE 

PRINT 975^ M/iMAXSQ 

975 FCRMaT 1/////1X,*THl NUNB-K OF INDEPENDENT OBSERVATIONS. ♦.I4.+. I 
IS GREATtR than THt NUMBER OF FUNCTIONS. *.IA.+. CONTAINED IN THE S 
2PECIFIFU MASTFk MODE L .♦ /I X . ♦ T HER EFOk t. THE UPPER FULL TRIANGLE OF 



SUBRDUTIKE SDBS2 7^/7^ QPT«1 


FTN 4*7+<»fl5 eO/Oi/23. 19.03.2< 


3THE MASTER MODEL WILL EL READ DIRECTLY FROM TAP- AND UStD IN THE 
_ <tFOLLDWING CALCULAT IONS . » ) 

300 CONTINUE 

130 INDEX ID— 999 

RETURN 

END 


RD HR. SEVERITY DETAILS 


D1A6F0SIS OF PROBLEM 


ro 

ro 



SYMBOLIC 

REFERENCE MAP 

(R-1) 





TRY 

POINTS 







3 

SUBS2 







ARIABLES SN 

type 

RELOCATION 






1 

INTEGER 


112 

IMAX 

integer 


113 

IMAXSQ 

INTEGER 


0 

INDEX 

INTEGER 

ARRAY 

141 

IROW 

INTEGER 


136 

IROWS 

iNTtGER 


125 

ISS 

INTEGER 


123 

ITEST 

INTEGER 


127 

ITS 

integer 


134 

K 

INTEGER 


135 

KStC 

INTEGER 


144 

K5TEP 

INTEGER 


140 

KTES 

INI EGER 


133 

KZON 

INTEGER 


137 

LEFT 

INTEGER 


0 

M 

INTtGER 


117 

MLIFF 

Integer 


C 

MiJRD 

INTEGER 


120 

MSEC 

integer 


122 

MTFS 

INTEGER 


121 

MZui't 

INTEGER 


115 

NISSEC 

INTEGER 


lib 

NISZON 

INTEGER 


ill 

NMSEC 

IMTEGEF. 


143 

NUMHTk 

INTEGER 


142 

NDMTP 

INTtGER 


126 

NSS 

INTEGER 


130 

Nssrc 

iNTcGtR 



F,P, 


F.P. 

F.P, 



5LilKDUTlH£ 5L6S3 


.74/74 OPT-1 


FTN 4*7+483 


eC/01/23. 19*03*20 


1 


5 


ID 


15 


20 

ro 

CO 


25 


30 


SUBRDUllNF SUBS3 (IND£X#NFUNC,ND£(j>«0R0 ) 

CQMF.fcNT — PARTICULAR CDtFFIClfcMTi OF INTfcREST ARE SPECIFIED AS CARD 
C INPUT TO BE READ FROM THIS SUBROUTINE* THE NUMBER OF 

C .FUNCTIONS^ NFUNCj TO BE HEAD IS A FORMAL PARAMETER OF TH13 

C SUBROUTINE AND DEFINES THE NUMBER OF ELEMENTS CONTAINED 

C IN THE ARRAY INDEX. 

DIMENSION INDtX(NFUNC) 

DO 100 i«l^NFUNC 
READ 15 * 915) LORO# LDE LEO 
COMMENT — LORD — ORDER OF FUNCTION* 

C LDtG — DEGREt Or FUNCTION* 

C LED — »0# FUNCTION IS EVEN. 

C =1, FUNCTION IS ODD. 

IF (LORD.EQ.O) GO TO 25 
IF (LDRD.cQ.LDEG) GD TO 50 

K«l NDEG+1 ) +2*MQRD+J.t0+L ORO+ (2+i.JEG-LORD-l )-l 
GO TU 75 
25 COFTINUE 
K-LDEG+1 
GD TD 75 
50 CONTINUE 

K»{KUEG+l)+2+LuRD-l+LED 
75 continue 

PRINT 925# i#LDRD#LDEG# L£D#K 
INLEX(1)»K 
ICO CONTINUE 
RETURN 

915 FORMAT (12# 12# ID 

925 FORMAT ( IX# +FUNCTION NUNbER +#14#* HAS BEEN SPECIFIED AS ORDER- *» 
112#*# CEGRtf *#I2#*# AND LED- *#I1#*. K- *»T4) 

ENL 


SYMBOLIC KEFtRENCE MaP (R»1) 

NTRY POINTS 

_ 3 sues3 



SLiikLLTIKL PRISYMI . 


74/74 


DPT-1 


PTU 4*7+465 


8C/01/23* 19*03.20 


. i 


5 


J.O- 


15 


20 


ro- 

■p> 


25 


30 


35 


40 


SUfckGUTlNE PRTSYMl ( S #NVV# NFUNC* ICODE ) 

COKHSNT — S IS A VECTOR WHICH CCNTAlNSTHt UPPER FULL TRIANGLE OF 
C SOME SYMMETRIC MATRIX 2* S IS PACKED AS FOLLOWS# 

C K»0 

C 00 20 J«l#tNFUNC 

C DO 10 I-1»J 

C K«K+1 

C S(KJ«Z(I»J) 

C 1C CONTINUE 

C 20 CONTINUE 

C WHERE NFUNC IS THF OrtDER OF THE MATRIX Z* 

C 

C PRTSYMl WILL EITHER# 

C A. PRINT THE VECTOR S AS THE UPPER FULL TRIANGLE OF Z# 

C ICCDE"o# Ok# 

C o» PRINT THE XQRRrLATION FORM OF ZCSAME fORMAT OIS OiflOVXJ — 

C FOR ICUDE-1, 

C 

C NFUNC AS DEFiNcD IN GLSRAN2 IS THE NUMBER OF COLUMNS IN THE 

C Z-MATRIX* 

C 

C MCDL ■ NUMBER Or (OLUKNS/LINE OF PRINT AND JIUST 3E IN 

C AGREEMENT WITH FORMAT STATEMENT 900 AND 950. 

C 

C THE ARRAY A DlMENSlUNED AS A(MCDL) IS USED FOR PRINTING 

C SO THAT S WILL NOT BE DcSTkOYEC WHEN ICDDE-1. 

C 

C Z iS DIVIDED INTO IR SECTIONS FOR PRINTING. 

C 

C JIDX - J INDEX OF FIRST COLUMN 

C FPF A PARTICULAR ScCflCN# KSEC. 

C JNDX - J index of last COLUMN 

C FOR A PARTICULAR SECTiON# KSEC. 

C JNDX ALSO « I INDFX OF LAiT ROW 

C FOR A FARTiCULAR SECTION# KSEC. 

C KSEC « number of SECTION BEING PRINTED. 

EUEI-SIDN S{NVV)#A(12) 

MCDL -12 
JNDX-C 

IRATIO-nFUNC/MCuL 

iR«IRATIO+i 

IDI FF«NFUNC-IRATI0+MC01 



SUbKDUTlNE 


45 


50 


55 


60 


65 


70 


75 


PRTSYMl 74/74 DPT-1 


FTN 4.7+485 80/01/23. 191.03^20 


If (IDIFf.N£.0) GO TO 25 

IP»IRAT10 

ICIFF-MCOL 

25 CONTINUE , 

DU 100 KS£C-1jIR 
JUX-JNtX + 1 
jNCX*JNtiX+NCJL 

IF (KSfcC.EO.IR) JNDX«( JNOX-MCOL )+IDIFF 

JEND«JNDX-(KStC-l)+MCOL 

PRINT 95(/j (J^J-JIDX,JNDX) 

Oil ICO I-l/JNDX 
DU 30 J-J1DX>JNDX 
J1«J-(KSEC-1)+MC0L 
A( JD-O. 

30 CONTINUE 

If (i.GT.JiUX) JIDX-JIDX+I 

DO 50 J»JltX#JNDX 

Jl»J“(KSEC-l)*hCOL 

K»( JX‘{ J-1) )/2+I 

If (ICuDE.LQ.O) GO TJ 35 

Kl«=(l+ (I+l) ) /2 

K2- { J* (J+1 J )/2 

A(ol)-5(K)/50RT (S(K1)*S (K2) ) 

60 TO 50 
25 CCKTJNUE 
A( J1)-S{K) 

50 CONTINUE 

PRINT 900. 1. (A( Jll #J1«I.J£ND) 
lUO CONTINUE 
RETURN 

900 FCPNAT (T4. 13.12(810.3)) 

950 fORNAT ( / /Til. 12 1 1 3. 7X ) / ) 

END 


SYNbOLlC RtFENENCF MaP (R-1) 



Appendix H - FOURIER SERIES REPRESENTATION OF A DISCRETE DATA SET 
The Fourier series approximation takes the form 

q 

g(x) = A + z [A.cos(£x) + B sin(£x)], (H-1) 

0 £ £ 

where g(x) is periodic over 2ir and q i Q. If now f(x) is a function for 
2Q + 1 discrete equally spaced values of x over the same period as g(x) above, 
a set of Fourier coefficients may be found that satisfies equation (H-1) by 
using the least-squares criterion. Let be the error associated with the 
rth value of x, then 

= f(Xp) - g(x^). 


and 

= Cf(x^) - Aq - E (Aj^cos(£x^) + B^sin(£x^))]^. (H-2) 

£— 1 

The least-squares technique as discussed earlier in this report requires 
that 


E[f(x„) - A- - I (A cos(£x J + B sin(£X„))]^ 


r' ”0 ", ' '£ 

r £=1 ^ 


r' £ 


be minimized. This leads to 


1 ^ 

A = z f(xp) > 

° r=-Q+l 


(H-3a) 


1 ^ 

An = jY E f(x„)cos(£X ), 

r=-Q+l 

for £ 7 < 0, Q, 


(H-3b) 


1 Q 

Afi "" m ^ f(x„)cos(Qx ), 

q r=-Q+l 


(H-3c) 
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and 


(H-3d) 
(H-4a) 

2 1 

= Q Hq + H-j cos(.£x-j) + H 2 cosCiiXg) + . . . + 

Hq_^ cos(£xq_-|) +^Hq cos(£Xq)], (H-4b) 

Aq = ^ [^Hq - + H 2 - . . . + (-1)^"'' Hq_^. (H-4c) 

and 

®£ ' f ■'■ ^2 sin(£X 2 ) + . . . + Gq_^ sin(£Xq_-|)], (H-4d) 

where 

H{x) = I [f(x) + f(-x)] (H-5a) 

and 

6(x) =l[f(x) - f(-x)]. (H-5b) 

Equations (H-4) and (H-5) then give the required Fourier coefficients which 
will satisfy equation (H-1). 

Since f(x) is periodic over 2 ti, 

f(ir) = f(-Tr) 

so that there are 2Q independent pieces of data. Then, as the objective of 
this Fourier series representation is data interpolation, rather than say 
"smoothing", the q in equation (H-1) takes on the value Q so that all 2Q 
possible terms are used. Equation (H-1) may now be written as 

Q 

f(x) = A^ + I [A^cos(£x) + B^sin(£x)]. (H-6) 

1 


®£ ^ f(x )sin(£x ). ^** 

^ ^ r=-Q+l ^ 

Equations (H-3) may be rewritten as 

^ V- 
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APPENDIX I - THE 0ZSTAT2 PROGRAM 


This appendix contains a listing of the OZSTAT2 program and its subroutines. 
The 0ZSTAT2 program has three basic capabilities: 

a. Data Grouping 

b. Statistical Analysis 

c. Computer Graphics 

The program is designed to read the BUV data formatted as described in 
Appendix B and to group the data into a global grid system. Based on 
this grid system means and variances are calculated for individual grid 
blocks, latitudinal zones, and for that part of the grid system that 
contains data. These calculations are described in more detail in section 
3. 

Graphics capabilities included in the 0ZSTAT2 program provide for each case 
a plot of the zonal means with ±1 sigma error bars, a scatter diagram of 
the ozone distribution as a function of latitude, and histograms of the 
data sampling distribution as a function of latitude or longitude. 






PROGRAM 

0ZSTAT2 74/7A 0PT»1 FTN 4,7+485 

80/01/23. 

14.08 

1 


PROGRAM 0ZSTAT2 { INPUT# 0UTPUT»TAPE5-INPUT»TAPE6«0UTPUT#TAPE10# 

□ ZSTAT2 

1 



1TAPE11#TAPE12#TAPE13#TAPE14#TAPE15#TAPE1#TAPE2) 

0ZSTAT2 

2 



COMMON/DO/DATE 

□ ZS TAT2 

3 



DIMENSION A(4)#K(5)#ALAT(100)#A0Z(100)#DATE(2) 

0ZSTAJ2 

. 4 

5 


DIMENSION SUMXSQ(36,24>#KK(36#24)#GP(36#24)#SUMX(36#24) 

QZSTAT2 

5 



DIMENSION KC0UNT(36#24»#SUMGPSQ(36#24)# VARGP(36#24) 

0ZSTAT2 

6 



DIMENSION ILAT(36)#RLAT(38)#RKK(36#24)#RLATBND(38)#NAHE(36) 

□ZSTAT2 

7 



DIMENSION RLNGBK(26),RK(26) 

0ZSTAT2 

8 



DIMENSION VARLATB(36)#AVGLATB(36)#RAT(36)»STDEVP(36),STDEVN{36) 

0ZSTAT2 

9 

10 


DIMENSION XDATA(6) 

□CT79 

1 



DIMENSION SUMT(36#24)#SUMTSQ(36,24)#SUMLT(36#24)#SUMLG(36#24) 

0CT79 

2 



DIMENSION SUMLTS0(36»24)#SUMLGSQ(36>24) 

0CT79 

3 



DATA NAME/3HI-1»3HI-2#3HI-3#3HI-4#3HI-5#3HI-6#3HI-7#3HI»8#3HI-9j4H0ZSTAT2 

12 



1I-10,4HI-11,4HI»12#4HI»13#4HI«14,4HI«15#4HI»16#4HI»17#4HI«18/ 

0ZSTAT2 

13 

15 

C 

SET UP PLOT VECTOR FILE 

0ZSTAT2 

14 


C 

SUBROUTINE PSEUDO (FN)# FN « FILENAME# CAN BE FOUND IN SECTION 

OZSTAIZ 

15 


C 

1.4.1 OF THE GRAPHICS MANUAL 

DZSTAT2 

16 


C 

IF LEROY IS NOT SPECIFIED# LIQUID INK PEN# BALL POINT PEN 

OZSTATZ 

17 


C 

IS AUTOMATICALLY CALLED# IF REQUIRED# BY DEFAULT. 

0ZSTAT2 

18 

20 

C 

LEROY IS ONLY USED FOR THE CALCOMP POSTPROCESSOR. 

0ZSTAT2 

19 


C 

FIRST FRAME MUST CONTAIN AT LEAST FIVE PLOT VECTORS WHEN USING 

0ZSTAT2 

20 

ro 

C 

CALCOMP. 

0ZSTAT2 

21 

UD 


CALL PSEUD0(6LMYSAV1) 

0ZSTAT2 

22 



CALL LEROY 

□ZSTAT2 

23 

25 


DO 10 I-l#6 

QZSTAT2 

24 



10 CALL CALPLT (0.#0.#-3) 

0ZSTAT2 

25 


C 

♦ ♦♦♦♦♦♦ 

0ZSTAT2 

26 


C 

INITIALIZE PARAMETERS FOR SUBROUTINE SELECT OPTIONS 

0ZSTAT2 

27 


c 

IF JSELECT -0 DO GRID BLOCKS ONLY 

0ZSTAT2 

28 

30 

c 

IF JSELECT -1 DO GRID POINTS ONLY 

0ZSTAT2 

29 


c 

IF JSELECT -2 DO BOTH 

0ZSTAT2 

30 



JSELECT-1 

0ZSTAT2 

31 



JSELECT-2 

0ZSTAT2 

32 



JSELECT-0 

0ZSTAT2 

33 

35 

c 

MSELECT -0# SKIPS BOTH PLOT ROUTINES 

0ZSTAT2 

34 


c 

MSELECT-1# CALL AVARPLT ONLY 

□ZSTAT2 

35 


c 

MSELECT-2# CALL HISTPLT ONLY 

DZSTAT2 

36 


c 

MSELECT -3# CALLS BOTH 

0ZSTAT2 

37 



MSELECT «0 

DZSTAT2 

36 

^0 


MSELECT -1 

0ZSTAT2 

39 



MSELECT -2 

0ZSTAT2 

40 



MSELECT -3 

0ZSTAT2 

41 




PROGRAM QZSTAT2 7A/7A DPT-1 FTN A.7+A85 

80/01/23. 

14.08.48 


C 

♦ ♦♦♦♦♦♦ 

□ZSTAT2 

42 


C 

INITIALIZE PARAMETERS FDR SELECTING GRID BLOCK SIZE 

0ZSTAT2 

43 


C 

ISIZE IS THE LATITUDE DIMENSION IN DEGREES 

0ZSTAT2 

44 


C 

JSIZE IS THE LONGITUDE DIMENSION IN DEGREES 

0ZSTAT2 

45 



ISIZE-5 

0ZSTAT2 

46 



JSIZE-15 

0ZSTAT2 

47 


C 

NLAT IS THE NUMBER OF ISIZE LATITUDE ZONES 

0ZSTAT2 

48 

50 

C 

NLONG IS THE NUMBER OF JSIZE LONGITUDE BANDS 

0ZSTAT2 

49 


C 

DIMENSION STATEMENTS MUST BE ADJUSTED FOR EACH RUN ACCORDING 

0ZSTAT2 

50 


C 

TO NLAT AND NLONG. 

DZSTAT2 

51 



NLAT-180/ISIZE 

0ZSTAT2 

52 



NL0NG-360/JSIZE 

0ZSTAT2 

53 

55 


PI»AC0S(-1. ) 

0CT79 

4 


C 

NCALC IS THE NUMBER OF TIME PERIODS OVER WHICH CALCULATIONS 

0ZSTAT2 

54 


C 

WILL BE EXECUTED. 

0ZSTAT2 

55 



NCALC-8 

0CT79 

5 



NCALC- 1 

0CT79 

6 

60 

C 


QZSTAT2 

57 



DO 200 L-l^NCALC 

0ZSTAT2 

58 



READ (5>225) KA#NDAY#DATE 

0ZSTAT2 

59 

Ca> 

c 

KA FROM TIME INTERVAL NCALC+1 MUST BE GREATER THAN NDAY 

0ZSTAT2 

60 

O 

c 

FROM TIME INTERVAL NCALC 

aZSIAT2 

61 

65 

c 

if 4< 4< 4< 1 4i 4c A * ♦ 4< ^ ^ ^ 4> 4< ♦ 

0ZSTAT2 

62 


c 

If 4c 4c If 4i 4c If * 4c 4c « 4c 4t If 4c 4> 4> « 4c 4> If 4> # If 4i 4nf # * « 4c « ^ If 4> 4> 4> # 4c « * 4nf # * * 4c ^ 4c « 4< « 4i 4c « « 4i 4i ft « 4< A « 

0ZSTAT2 

63 


c 

TOTAL OZONE DATA IS AVERAGED. MEAN AND VARIANCE ARE PUT INTO 

0ZSTAT2 

64 


c 

PARTICULAR ELEMENTS OF AN NLAT X NLONG GRID SYSTEM. ♦4c4c 

0ZSTAT2 

65 


c 

4 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c # 4c 4c 4c 4< 4c 4c if 4c If 4c 4c » 4c 4> 4c 4c 4c 4c 4c 4c 4c If 4c If 4c 4c 4> # 4c ♦ 4c 

0ZSTAT2 

66 

70 

c 

IN THE FOLLOWING STATEMENTS# CERTAIN PARAMETERS ARE INITIALIZED 

0ZSTAT2 

67 


c 

IJ IS USED AT STAEMENT 57 

0ZSTAT2 

68 



IJ— 1 

0ZSTAT2 

69 


c 

lEOF IS USED AT STATEMENT 21. IEOF-1 INDICATES THAT THE 

0ZSTAT2 

70 


c 

END OF FIL-E HAS BEEN REACHED. 

0ZSTAT2 

71 

75 


I^OF-0 

DZSTAT2 

72 



IDAY-K4 

0ZSTAT2 

73 



IISUM-0 

DZSTAT2 

74 



GSUM-0. 

□ZSTAT2 

75 



DO 15 I»1,NLAT 

0ZSTAT2 

76 

80 


DO 15 J-1#NL0NG 

0ZSTAT2 

77 



KK(I,J)«0 

DZSTAT2 

79 



SUMXSQd# J)-0. 

0ZSTAT2 

80 



SUMX(I#J)«0. 

0ZSTAT2 

81 



SUMTU, J)-0. 

0CT79 

7 



» 


PROGRAM QZSTAT2 7A/7A OPT-1 FTN A.7+A85 80/01/23. 14,08.48 


85 



SUMTSQ(I#J)-0. 

0CT79 

8 




SUMLTd, J)-0. 

0CT79 

9 




SUHLG(1> J)>0. 

0CT79 

10 




SUMLTSQ(IjJ)-0. 

0CT79 

11 




SUHLGSQ<I> J)»0. 

OCT 79 

12 

90 


15 

CONTINUE 

0ZSTAT2 

82 


C 



0ZSTAT2 

83 




IF( (MSELECT.EO.O).OR. (MSELECT.EQ.2) ) GO TO 18 

0ZSTAT2 

84 


C 


CALL AVARPLT TO CONSTRUCT AXES AND LABELS FOR SCATTER DIAGRAM, 

0ZSTAT2 

8 5 


C 


MEAN CURVE AND STANDARD DEVIATION PLOT. DEFINE M AS ANYTHING. 

0ZSTAT2 

86 

95 



M«100 

0ZSTAT2 

87 




CALL AVARPLT ( VARLATB,AVGLATB,ALAT, AOZ,RAT,STDEVP,STDEVN,NLAT) 

0ZSTAT2 

88 


C 



0ZSTAT2 

89 



18 

CONTINUE 

0ZSTAT2 

90 




DO 35 M-1,100 

0ZSTAT2 

91 

100 


20 

READd) (XDATA( J)#J-1#6) 

0CT79 . 

13 




IF (EOF(D) 21,22 

0ZSTAT2 

93 



21 

IEOF-1 

0ZSTAT2 

94 




GO TO 23 

0ZSTAT2 

95 



22 

CONTINUE 

0ZSTAT2 

96 

105 



K(4)»XDATA(2) 

0CT79 

14 




K(4).(XDATA(1)-1970. )*365.+K(4) 

Q£T79^ 

15 

CO 



IF (XDATA(1).EQ.1973..0R.XDATA(1) .EQ.1974..0R.XDATA(1).EQ.1975.. 

0CT79 

16 



IDR.XDATA(l) .EQ.1976.) K(4)-K(4)+l 

0CT79 

17 




IF (XDATAd ) .E0.1977. ) K(4)»K(4) + 2 

OCT 79 

18 

110 



RK5»XDATAC3)/3600. 

0CT79 

19 




K{5)-RK5 

0CT79 

20 




A(1)«(RK5-K{ 5) )+60. 

OCT 7.9 . 

21 




A(2)«X0ATA(4) 

0CT79 

22 




A(3)«XDATA(5) 

0CT79 

23 

115 



A(4)-ABS(XDATA(6) ) 

0CT79 

24 




IF (A(4).EQ.999,,DR.A(4).EQ.77.) GO TO 20 

0CT79 

25 


C 


KC4) IS THE DAY NUMBER ♦♦♦♦♦♦♦♦♦♦*♦♦♦♦♦+♦♦♦♦ 

025 TAT 2 

97 




IF (A(4) ,LT. 0.200. OR. A(4).GT. 0,65) PRINT 905, XOATA 

0CT79 

26 



905 

FORMAT ( ♦0*,T10,6(E15.8, 5X)/ ) 

0CT79 

27 

120 



IF (K(4) .LT.IDAY) GO TO 20 

02STAT2 

98 




IF (K(4) .LE.NDAY) GO TO 30 

0ZSTAT2 

99 




BACKSPACE 1 

0ZSTAT2 

100 



23 

IF (M.EO.l) GO TO 25 

0ZSTAT2 

101 




IF( (MSELECT.EO.O). OR. (MSELECT.EQ.2) ) GO TO 25 

0ZSTAT2 

102 

125 



CALL SCAT ( VAPLATB, AVGLATB,ALAT,A0Z,M-1,RAT,STDEVP,STDEVN,NLAT) 

0ZSTAT2 

103 



25 

IF (lEOF.EQ.l) GO TO 75 

0ZSTAT2 

104 



PROGRAM 02STAT2 


74/74 OPT-1 


FTN 4.7+485 


80/01/23. 14.08,48 





GO TO 50 

0ZSTAT2 

105 


C 


A(2) IS THE LATITUDE 

0ZSTAT2 

106 



30 

ALAT(H)>A(2) 

0ZSTAT2 

107 

130 



A0Z(H)«A<4) 

0ZSTAT2 

108 




K4»K{4) 

□ZSTAT2 

109 


C 


0f^^4*******^***tlf********i^**i^0>lt***********4f**if*****ti************* 

0ZSTAT2 

110 


C 


DATA RECORDS FOR WHICH THE LATITUDE - 0 DEGREES WILL BE ASSIGNED 

QZSTAT2 

111 


C 


TO THE NORTHERN HEMISPHERE. HOWEVER# IT IS BELIEVED THAT NO 

0ZSTAT2 

112 

135 

C 


SUCH RECORDS OCCUR BETWEEN DAYS 101 AND 465. 

0ZSTAT2 

113 


C 


* 4< * « « * « 4^ 4> « 4> 4 4 4 4 4< >4’ 4> 4t 4 « 4< 4 4 4 4 >t> 4 4> 4 4 4 4 m 4< 4n<> 4 # 4 4 

0ZSTAT2 

11^ 




LAT-ABS(A(2)) 

DZSTAT2 

115 




I-LAT/ISIZE+1 

0ZSTAT2 

116 




IF (A(2).LT.O) I-I+NLAT/2 

0ZSTAT2 

117 

140 

C 


AC3) IS THE LONGITUDE 

0ZSTAT2 

116 




L0NG>A(3) 

□ZSTAT2 

119 




J-LONG/JSIZE+1 

OZSTATZ 

120 




KK(I#J)»KK(I#J)+1 

0ZSTAT2 

129 




XSQ-A<4)4A(4) 

0ZSTAT2 

130 

145 



SUHXSO(I»J)-SUMXSO(I# J)+XSO 

0ZSTAT2 

131 




SUMXd, J)»SUMX(I# J)+A(4) 

0ZSTAT2 

132 

to 



TM«K(4)+XDATA(3)/86400. 

0CT79 

28 

ro 



TMSQ-TH4TM 

0CT79 

29 




SUHTSOd# J)«SUMTSQ{I# J)+TMSQ 

0CT79 

30 

150 



SUMT(I#J )-SUMT(I»J )+TM 

0CT79 

31 




SUMLTCI# J)-SUMLT(I#J )+A(2) 

□ CT79 

32 




SUHLGd# J)«SUMLG(I> J)+A(3) 

0CT79 

33 




SUMLTSQd# J )-SUMLTSQ(I» J )+A(2)4A(2) 

0CT79 

34 




SUMLGSQd# J)-SUMLGSQ(l# J)+A(3)4A<3) 

QCT79 

35 

155 


35 

CONTINUE 

0ZSTAT2 

133 




IF( (MSELECT.EO.O).OR. (MSELECT.EQ.2) ) GO TO 18 

□ZSTAT2 

134 




CALL SCAT {VAPLATB»AVGLAT8#ALAT#A0Z#100#RAT»STDEVP#STDEVN#NLAT) 

0ZSTAT2 

135 




GO TO 18 

0ZSTAT2 

136 


C 


4444444444444444444444444444444444444444444444444444444444444444 

0ZSTAT2 

137 

160 

C 


444444444444444444444444444444444 

0ZSTAT2 

138 


c 


BEGIN STATISTICS CALCULATIONS FOR GRID SYSTEMS. 

0ZSTAT2 

139 



50 

IF (JSELECT.EO.l) GO TO 97 

0ZSTAT2 

140 




PRINT 115# NLAT#ISIZE# JSIZE#NLAT/2#NLAT/2+l#NLAT#NLQNG 

□ZSTAT2 

141 




PRINT 116# DATE 

QZSTAT2 

142 

165 



PRINT 117# L 

0CT79 

36 




GVAR-0. 

OCT 79 

37 




KBLK-O 

0CT79 

38 




DO 65 I-1#NLAT 

0ZSTAT2 

143 




PROGRAM 

0ZSTAT2 7A/7^ DPT-1 FTN 4.7+<»85 

80/01/23. 

14.08.48 




SUM-0. 

0ZSTAT2 

144 

170 



ISUM-0 

0ZSTAT2 

145 




SSQLATB-0. 

0ZSTAT2 

146 




TSUM-0. 

OCT 7.9 

39 




TSQLATB-0. 

QCT79 

40 




PRINT 120 

0ZSTAT2 

147 

175 

C 



02STAT2 

148 




DO 55 J-1,NL0NG 

0ZSTAT2 

149 




IF (KK(I^ J}.EQ.O) GO TO 51 

DZSTAT2 

150 




AX-SUMX(I#J}/KK{I> J) 

OZSIAT-2 

058 




AT-SUMT(I,J)/KK(IiJ) 

0CT79 

41 

180 



AVLAT-SUMLTd# J )/KK(I>J ) 

OCT 79 

42 




AVLDNG-SUMLGd^ J)/KK(I, J) 

OCT 79 

43 




IF <KK(I, J) .EO.l) GO TO 52 

DZSTAT2 

159 




VARX-(SUMXSQn»J)~KK(I» J )*AX*AX)/ (KK(I*J)-1. ) 

0ZSTAT2 

160 




VART-(SUMTSQ(I# J)-KK(I, J)*AT*AT)/CKK(I> J}-1.) 

flCT.7.9 

44 

185 



GO TO 53 

Q2STAT2 

161 



51 

AX-0. 

□ZSTAT2 

162 




AT-0. 

OCT 79 

45 




AVLAT-0. 

QCT79 

46 




AVLONG-0. 

0CT79 

47 

190 


52 

VARX-0. 

02SIilT2 

163 




VART-0. 

0CT79 

48 

LO 


53 

CONTINUE 

0ZSTAT2 

164 

CO 

C 


SUM IS ACCUMULATIVE OZONE CONTENT /LATITUDE BAND 

02STAT2 

165 




SUM-SUM+SUMX(I# J) 

0ZSTAT2 

166 

195 



TSUM-TSUM+SUMT(I,J ) 

0CT79 

49 


C 


GSUM IS THE GLOBAL ACCUMULATIVE OZONE LAYER THICKNESS 

Q2STAI2 

167 




GSUM-GSUM+SUMX(I» J) 

OZS TAT 2 

168 


C 


ISUM IS THE NUMBER OF DATA POINTS AVERAGED/LATITUDE BAND 

0ZSTAT2 

169 




ISUM-ISUM+KK(I,J) 

0ZSTAT2 

170 

200 



SSOLATB-SSQLATB + SUMXSQdi J) 

DZSTAT2 

171 




TSQLATB-TSQLATB+SUMTSQ(1#J) 

0CT79 

50 




PRINT 100# AX#VARX»KK(I# J),I#J,AVLAT#AVLONG#AT»VART 

0CT79 

51 




IF (KK (I# J ) .EO.O) GO TO 55 

0CT79 

52 




KBLK-KBLK+1 

DCT79 

53 

205 



WRITE (10) (90,-AVLAT)*PI/180.#AVL0NG*PI/180,# AX^IOOO. 

0CT79 

54 


C 



0CT79 

55 


C 


♦ ♦♦♦♦♦♦♦♦ 

0CT79 

56 


C 


XLATd# J)-AVLAT 

0CT79 

57 


C 


XLONGd# J)-AVLONG 

0CT79 

58 

210 

C 


XOZd» J)-AX 

0CT79 

59 



215 


220 


225 


230 


OJ 

235 


2A0 


2A5 


250 


PROGRAM 0ZSTAT2 74/74 OPT-1 



FTN 4. 

7+485 

80/01/23. 

14.08.48 

C 


♦ ♦ ♦ 4 ]|c 

4 

4 

4 


0CT79 

60 

C 


♦ ♦ ♦ ♦ ^ 

4 

4 

4 


0CT79 

61 

C 


GRID BLOCK DATA FDR NCALC WEEKS IS 

STORED 

IN 

ZDATA. 


0ZSTAT2 

174 

C 


NUMBER OF DATA POINTS/GRID BLOCK FOR NCALC 

WEEKS 


0ZSTAT2 

175 

C 


IS STORED IN NDATA. 





0ZSTAT2 

176 


55 

CONTINUE 





0ZSTAT2 

179 

C 







0ZSTAT2 

180 

C 


IISUM IS THE TOTAL NUMBER OF OZONE 

RECORDS 

USED 


0ZSTAT2 

181 

C 


IN THESE CALCULATIONS. 





DZSTAT2 

182 



IISUH-IISUM+ISUM 





0ZSTAT2 

183 

C 


<M> * « 4> ♦ 4< 4> 4< 4> >•< « 4< « 4: Hi ♦ * « # 4> 4i 4> « « 4> « ♦ >l> « 4c 4c 



0ZSTAT2 

184 

c 


ILAT(I), I-l,NLATf IS THE NUMBER OF 

DATA POINTS/LATITUDE BAND, 

0ZSTAT2 

185 

c 


SUCH THAT I-1»NLAT CORRESPONDS TO LATITUDE 

BANDS 

FROM SOUTH 

0ZSTAT2 

186 

c 


TO NORTH. ILAT IS AN INPUT PARAMETER OF 

HISTPLT. 


0ZSTAT2 

187 



IF (I.GE.NLAT/2+1) GO TO 57 





0ZSTAT2 

188 



ILAT(I+NLAT/2)-ISUM 





0ZSTAT2 

189 



GO TO 58 





0ZSTAT2 

190 

c 


IJ IS INITIALIZED AS -1. 





0ZSTAT2 

191 


57 

IJ-IJ+2 





0ZSTAT2 

192 



ILAT(I-I J)-ISUM 





0ZSTAT2 

193 


58 

CONTINUE 





aZSTAT2 

194 

c 


^ii4‘*********’t‘************************m**t*** 



0ZSTAT2 

195 



IF (ISUM.EQ.O) GO TO 60 





0ZSTAT2 

196 



AX-SUM/ISUM 





0ZSTAT2 

197 



AT-TSUM/ISUM 





0CT79 

62 



GVAR-GVAR+SSQLATB 





0CT79 

63 



IF (ISUM.EQ.l) GO TO 61 





0ZSTAT2 

198 



VARX-(SSQLATB-AX*AX*ISUM)/(ISUM-1.) 





_Q2STAT2 

199 



VART«(TSQLATB-AT+AT+ISUM)/(ISUM-1.) 





0CT79 

64 



GO TO 64 





0ZSTAT2 

200 


60 

AX»0. 





0ZSTAT2 

201 



AT-0. 





0CT79 

65 


61 

VARX-0. 





0ZSTAT2 

202 



VART-0. 





0CT79 

66 


64 

VARLATB<I)«VARX 





0ZSTAT2 

203 



AVGLATB(I)-AX 





0ZSTAT2 

204 

c 


4 4 4 4 4 4 

4 

4 

4 

4 

DCT79 

67 

c 


WRITE (10^910) AX41000.,SQRT{VARX41000000. 

) 



0CT79 

68 

c 

910 

FORMAT (1X,2(F7.3#2X ) ) 





0CT79 

69 

c 


4 4 4 4 4 4 

4 

4 

4 

4 

0CT79 

70 


65 

PRINT 110, I,AX,ISUM, VARX,AT,VART 





0CT79 

71 

c 


4 4 4 4 4 4 

4 

4 

4 


0CT79 

72 



PROGRAM 0ZSTAT2 


74/74 OPT-1 


FTN 4.7+485 


80/01/23. 14.0^.48 



C 


♦ ♦♦♦♦♦♦♦♦ 

0CT79 

73 


c 


WRITE(IO) KK 

0CT79 

74 

255 

c 


WRITE(IO) XLAT 

0CT79 

75 


c 


WRITE(IO) XLONG 

I3I:T79 

76 


c 


WRITE(IO) XOZ 

0CT79 

77 


c 



QCT79 

78 


c 


* .4> t « * * * 4> 

0CT79 

79 

260 

c 


# 4c « 4c « 4ntr « * 4i 4t 4c 4t 4^ 4> « 4> 4< 4c # 4> * # 4: 4> 4< 4< 4i « 

0ZSTAT2 

206 


c 



DZSTAT2 

207 




PRINT 125> IISUMjIDAYiK4 

QZSTAT2 

208 




GAV-GSUM/IISUM 

0ZSTAT2 

209 




PRINT 130#GAV 

JJZS TAT2 

210 

265 



GVAR«(GVAR-GAV4cGAV*IISUM)/(IISUM-l.) 

0CT79 

80 




PRINT 135, GVAR,SQRT(GVAR) 

0CT79 

81 



135 

FORMAT (1X,4«gL0BAL VARIANCE- ♦,E15.8/ 

0CT79 

62 




llX,4cSTAN0ARD DEVIATION- 4.,E15.8) 

ACT 79 

82 




PRINT 925, KBLK 

□CT79 

84 

270 


925 

FORMAT C4c0*,*A TOTAL OF ♦,I6,4c BLOCKS ARE FILLED.*) 

0CT79 

85 


c 


4c44c4<«4c<4444ti4*4i4c«4>4i44c*4c4c4c4c4c4‘4c4c4c4cc44c4c4>4c#4c4c4c4c4>«4c4>*4:4c4>*4i4c4>4c4‘****4c**4>*4c 

0ZSTAT2 

211 




IF (MSELECT.EO.O) GO TO 71 

0ZSTAT2 

212 




IF (MSELECT.EQ.2) GO TO 70 

0ZSTAT2 

213 




CALL ASTD ( VARL ATB, AVGLATB, ALAT, AOZ,M, RAT,STDE VP,STOEVN, NLAT ). 

□ZSJAI2 

.214 

275 



IF (MSELECT.EO.l) GO TO 71 

0ZSTAT2 

215 



70 

CALL HISTPLT ( ILAT,KK,NLAT,NLONG,RLAT,RKK,RLATBND,RLNGBK,RK,NAME, 

0ZSTAT2 

216 




lNLAT+2,NL0NG+2) 

0ZSTAT2 

217 



71 

CONTINUE 

0ZSTAT2 

218 


c 


4c 4 4c 4c 4c 4c 4c 4< 4 4> 4t 4c 4< 4c 4> 4c 4> 4> ^ « 4c 4c 4> 4c 4> 4< 4c « 4< * * 4c 4> « « 4> 4c 4i 4c « 4> * 

0ZSTAT2 

219 

280 



IF ( JSELECT.EQ.2) GO TO 97 

AZSTAT2 

220 




GO TO 200 

0ZSTAT2 

221 



75 

PRINT 145, K(4) 

0ZSTAT2 

222 




IF(K(4) .EQ.173) GO TO 50 

0CT79 

86 




PRINT 160 

0ZSTAT2 

224 

285 



STOP 1 

0ZSTAT2 

225 



97 

CALL GRID (SUMX,SUMXSQ,KK,GP,NLAT,NLONG,KCOUNT,SUMGPSQ,VARGP) 

0ZSTAT2 

226 


c 


4c44444i4c4c44c4c4c4c*44c4c444c4*4c4‘4**4c4c4c4***4c4c4c44c4c4c4* 

0ZSTAT2 

227 


c 


4c 4 4c 4c 4 4 4c 4 4 4c 4c 4c 4 ** 4c 4c 4 4c 4c 4 4c 4c 4 4c 4c 4c 4c 4c 4 4c 4> 4c 4c 4c 4c 4 4c 4 4c 4c 4c 4 

QZSTAT2 

228 


c 


44444 4444444 4444444444444*44444**444444444* 

0ZSTAT2 

229 

290 

c 


4444444444444444444444444444444444444444444 

QZSTAT2 

230 



200 

CONTINUE 

0ZSTAT2 

231 




CALL NFRAME 

0ZSTAT2 

232 




CALL CALPLT (0.,0.,999) 

0ZSTAT2 

233 




STOP 

0ZSTAT2 

237 



FTN 4.7+485 


PROGRAM DZSTAT2 74/74 OPT-1 


80/01/23. 14.06.48 


295 


300 


305 


310 


315 

OJ 

CTi 


320 


100 FORMAT (1X#2(E16.8>5X)» I3i 5X#2(I2#5X)»T70,F7.3>T88#F8.3#T105#F8.3,DCT79 
1T119#E15.8) DCT79 

110 FORMAT l*0+i*THE AVERAGE OZONE DENSITY FOR THE LATITUDE BAND*/1X,*QZSTAT2 
ICORRESPQNDING TO I- IS ♦»E16.8#* THIS CALCULATION IS BASED 0ZSTAT2 

20N ♦,I6,* RECORDS OF DATA*/1X,*THE VARIANCE OF THE MEAN IS *, E16 . 80ZS TAT2 
1/1X#*TIME AVERAGE- F 8 . 3» 10X#*TIME VARIANCE- */E15.8) DCT79 

115 FORMAT (♦1*>+IN THE FOLLOWING TABLES^ OZONE DENSITY HAS BEE0ZSTAT2 

IN AVERAGED BY GRID BL OCXS . ♦/IX# ♦THESE BLOCKS REPRESENT AN AREA ON 0ZSTAT2 
2 THE EARTH"S SURFACE THAT IS DEGREES LATITUDE BY DEG0ZSTAT2 

3REES LONGITUDE. ♦/IX, ♦latitude BANDS CORRESPONDING TO LATITUDE INDIQZSTAT2 
4CES 1 THROUGH *p12p* ARE IN THE NORTHERN HEHISPHERE«/1X»^WHILE LAT0ZSTAT2 
5ITUDE INDICES *pI2p* THROUGH ♦,I2/^ ARE IN THE SOUTHERN HEMISPHERE0ZSTAT2 
6. LONGITUDE INDICES (1-^jI2#^) RANGE FROM^/IX/^THE GREENWICH HERI0ZSTAT2 
7DIAN WESTWARD THROUGH 360 DEGREES.*) 0ZSTAT2 

116 FORMAT (1X#*THESE CALCULATIONS ARE FOR THE TIME PERIOD ♦#A10#A4) QZSTAT2 

117 FORMAT (11X#*AND ARE FOR CASE NUMBER ♦^IZ) 0CT79 

120 FORMAT ( *0*» T8# *ME AN*, T27, ♦ VARI ANCE*» T43> *K < J ) ♦# T53# ♦!♦/ T60/ ♦ J*/ OZS TAT2 

1T70,*MEAN LATITUDE*, T88,*MEAN LONGITUDE*, T105,*AVERA6E TIME*, 0CT79 


2T119,*VARIANCE*> 

125 format (*0*,I7,* 

1S*/1X,*THIS DATA 
2 INCLUSIVE.*) 

130 FORMAT (*0*,*THE 

IE PERIOD IS *,E15.8) 

145 FORMAT (*0*,*REACHED END OF FILE. 


0CT79 

RECORDS OF DATA ARE USED IN THE ABOVE CALCULATI0N0ZSTAT2 
INCLUDES RECORDS FROM DAYS *,I4,* THROUGH *,I4,*,0CT79 

0ZSTAT2 

GLOBAL OZONE LAYER THICKNESS AVERAGE FOR THIS TIM0ZSTAT2 

0ZSTAT2 

DAY NUMBER - *,I5) QCT79 


160 FORMAT 
225 FORMAT 
END 


( ///*0*,*REACHED 
(2I4,A10,A4) 


EOF PRIOR TO LAST DAY ON FILE*) 


0ZSTAT2 

0CT79 

QZSTAI2 


87 

88 

239 

240 

241 

89 

243 

244 

245 

246 

247 

248 

249 

250 

251 

90 

252 

91 

92 
254 

93 

256 

257 

258 

94 
260 

95 
262 


ID NR. SEVERITY DETAILS DIAGNOSIS OF PROBLEM 

13 I NAME DATA VARIABLE LIST EXCEEDS ITEM LIST, EXCESS VARIABLES NOT INITIALIZED. 


SYMBOLIC REFERENCE MAP (R-1) 



I 


SUBROUTINE GRID 74/74 OPT-1 FTN 4.7+485 80/01/23. 14.08.48 


1 


SUBROUTINE GRID ( S# SS/K# GP# NLAT# NLONG# KCOUNT# SUMGPSQ# VARGP ) 

GRID 

1 



COMMON/DD/DATE 

GRID 

2 



DIMENSION S(NLAT>NLONG)#K(NLAT,NLONG)>GP(NLAT>NLONG)#KCOUNT(NLAT» 

NGRID 

3 



1L0NG)jSS (NLAT>NLONG)»SUMGPSQ(NLAT>NLONG)>VARGP(NLAT#NLONG) 

GRID 

4 

5 


DIMENSION DATE (2)' 

GRID 

5 



M-NLAT/2-1 

GRID 

6 



Ml-NLAT/2+1 

GRID 

7 



M2-NLAT/2 

GRID 

6 



M3-NLAT/2+2 

GRID 

9 

10 


PRINT 105» NLAT#M2,M1^M3»NLAT»NL0NG 

GRID 

10 



PRINT 110,DATE 

GRID 

11 


C 

4< 1# 4c 4i 4: 4i 4c # 4 e # 4c % 4c 

GRID 

12 


C 

GRID POINTS IN THE NORTHERN HEMISPHERE ARE CALCULATED BELOW 

GRID 

13 



DO 25 I»1,M 

GRID 

14 

15 


DO 25 J«1,NL0NG 

GRID 

15 



L-J+1 

GRID 

16 



IF (L.EQ.NLONG+1) L-1 

GRID 

17 



N-I+1 

GRID 

18 



SUMGPSQ (I#J)»SS(I»J )+SS(I>L )+SS(Nj J)+SS(N,L ) 

GRID 

19 

20 


KCOUNTd# J)-K{I,J)+K(I>L)+K(N» J)+K(N>L) 

GRID 

20 



IF (KCOUNTCI# J).E0.0) GO TO 20 

GRID 

21 



GP(I, J)-(S(I»J)+S{I#L)+S<N»J)+S(N»L))/KCDUNT{I*J) 

GRID 

22 

CO 


IF (KCOUNT(I#J).EQ.l) GO TO 21 

GRID 

23 

“vl 


VARGP (I# J)-(SUMGPSQ( I# J )-KCOUNT ( I# J) ♦GP ( I# J )*GP (I» J ) ) / (KCOUNT (1# JGRI D 

24 

25 


1)-1. ) 

GRID 

25 



GO TO 25 

GRID 

26 



20 GP(I#J)»0, 

GRID 

27 



21 VARGP (I,J) »0. 

GRID 

28 



25 CONTINUE 

GRID 

29 

30 

C 

4c 4i 1 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c + 4c 4c 4c 4c 4c 4c 4c 4c + 4c 4c 4c 4c 4< 4c 4c ♦ 4c 

GRID 

30 


C 

GRID POINTS ALONG THE EQUATOR ARE CALCULATED BELOW 

GRID 

31 



DO 30 J-1#NLQNG 

GRID 

32 



L-J+1 

GRID 

33 



IF (L.EQ.NLONG+1) L-1 

GRID 

34 

35 


SUMGPSQ(M1>J )-SS(l,J )+SS(l,L)+SS(Ml#J )+SS (M1#L ) 

GRID 

35 



KCOUNT (M1,J)-K(1,J)+K(1#L)+K(M1, J)+K(M1,L) 

GRID 

36 



IF (KC0UNT(M1, J).EQ.O) GO TO 27 

GRID 

37 



GP(M1# J)-(S(1# J) + S(l> L) +S { Ml# J )+S(Ml» LI) /KCOUNT (M1>J) 

GRID 

38 



IF (KCOUNKMl, J),EQ.l ) GO TO 28 

GRID 

39 

40 


VARGP (Ml# J)-(SUMGPSQ(M1# J)-KC0UNT(H1# J)*GP(M1# J)*GP(M1» J) )/(KCOUNGRID 

40 



1T(M1#J)-1.) 

GRID 

41 



GO TO 29 

GRID 

42 



SUBROUTINE 

GRID 

74/74 OPT-1 FTN 4.7+485 

80/01/23. 

14.08.48 



27 

GP (Hl^J)-O^ 

GRID 

43 



28 

VARGP (M1>J)»0. 

GRID 

44 

^5 


29 

GP (M2,J)-0. 

GRID 

45 




VARGP (M2»J)»0. 

GRID 

46 




KCQUNT <M2iJ)«0 

GRID 

47 



30 

CONTINUE 

GRID 

46 


C 



GRID 

49 

50 

C 


GRID POINTS IN THE SOUTHERN HEMISPHERE ARE CALCULATED BELOW 

GRID 

50 




DO 35 I-M3#NLAT 

GRID 

51 




DO 35 J-1,NLQNG 

GRID 

52 




L-J + 1 

GRID 

53 




IF (L.EQ.NLONG+1) L-1 

GRID 

54 

55 



N«I-1 

GRID 

55 




SUMGPSO (I,J)»SS(I^ J)+SS<I»L)+SS(N,J)+SS(N»L) 

GRID 

56 




KCOUNTd# J)»K(N, J )+K (N> L ) +K ( 1, J } +K C I# L ) 

GRID 

57 




IF <KCOUNT(I»J).EQ.O) GO TO 33 

GRID 

58 




GP(I#J)«(S(I»J)+S(I#L) + S(N#J )+S(N»U)/KCOUNT(I>J) 

GRID 

59 

60 



IF (KCOUNT(I»J).EQ.l) GO TO 34 

GRID 

60 




VARGP (I» J)-CSUMGPSQ( I# J )-KCOUNT ( I» J) ♦GP C J )*GP (I» J ) ) / ( KCDUNTCl# JGRI D 

61 



1 )-l. ) 

GRID 

62 




60 TO 35 

GRID 

63 



33 

GP(I»J)-0. 

GRID 

64 

65 


34 

VARGP (I,J)-0. 

GRID 

65 

CD 


35 

CONTINUE 

GRID 

66 


c 



GRID 

67 


c 


KC IS A COUNTER FOR THE TOTAL NUMBER OF DATA POINTS USED IN 

THESE GRID 

68 


c 


CALCULATIONS 

GRID 

69 

70 

c 


OZDEN IS AN ACCUMULATIVE SUM OF TOTAL OZONE DENSITY# SUMMED 

IN A GRID 

70 


c 


PARTICULAR LATITUDE BAND. 

GRID 

71 


c 


KCl IS A COUNTER FOR THE NUMBER OF DATA POINTS USED IN THE 

GRID 

72 


c 


CALCULATIONS FOR A PARTICULAR LATITUDE BAND. 

GRID 

73 


c 


*«««**« 1(1* 

GRID 

74 

75 



KC-0 

GRID 

75 




DO 55 I-1»NLAT 

GRID 

76 




OZDEN-0. 

GRID 

77 




KCl-0 

GRID 

78 




PRINT 100 

GRID 

79 

80 



DO 50 J»1#NL0NG 

GRID 

80 




KC»KC+KCOUNT(I#J) 

GRID 

81 




KC1-KC1+KC0UNT(I#J ) 

GRID 

82 




OZDEN-OZDEN+GP(I» J )+KCOUNT(I#J) 

GRID 

83 



50 

PRINT 125# GP(I#J)#KCOUNT ( I# J ) > I# J# VARGP II#J) 

GRID 

84 


• # « . . .j. » 



SUBROUTINE GRID 


7^/7<t OPT-1 


FTN 4.7+<»85 


80/01/23. 14.08.48 


85 


90 


95 


100 


105 


AVDEN-QZDEN/KCl GRID 85 

55 PRINT 150, I,AVDEN,KC1 GRID 86 

PRINT 175,KC GRID 87 

100 FORMAT <*0+,T4,*MEAN OZONE DENSITY*, T27, ♦VARIANCE*, T43,*KC0UNT*,T5GRID 88 

13,*I*,T60,*J*) GRID 89 

105 FORMAT (*1*,*IN THE FOLLOWING *,I2,* TABLES, OZONE DENSITY HAS BEEGRID 90 

IN AVERAGED. BY GRID POINTS. LATITUDE INDICES 1 THROUGH *,I2 GRID 91 

2/lX,*C0RRESP0ND TO THE NORTHERN HEMISPHERE, LATITUDE INDEX *,I2,* GRID 92 

3CORRESPONDS TO THE EQUATOR, AND LATITUDE INDICES *,I2,* THROUGH *,GRID 93 

4I2/1X,*C0RRESP0ND TO THE SOUTHERN HEMISPHERE. LONGITUDE INDICES AGRID 94 

5RE SIMILAR TO THOSE ABOVE WITH J- *,I2,* CORRESPONDING TO LONGITUOGRID 95 

6E«0*/1X,*DEGREES OR THE MERIDIAN THROUGH GREENWICH*) GRID 96 

110 FORMAT (1X,*THESE CALCULATIONS ARE FOR THE TIME PERIOD *,A10,A4) GRID 97 

125 FORMAT ( IX, E16 . 8, T44, 14, T53, 12, T60, 12, T25, El 6 . 8 ) GRID 98 

150 FORMAT (*0*,*THE AVERAGE OZONE DENSITY CORRESPONDING TO LATITUDE IGRID 99 

INDEX *,I2,* IS ♦, E16.8,*. THIS MEAN IS BASED ON *,I6,* DATA POINTGRID 100 

2S.*) GRID 101 

175 FORMAT (*0*,I7,* DATA RECORDS ARE USED IN THE ABOVE CALCULATIONS .*GRID 102 

1) GRID 103 

RETURN GRID 104 

END GRID 105 


CjO 

LO 


SYMBOLIC REFERENCE MAP (R-1) 

TRY POINTS 
3 GRID 


RIABLES 

SN TYPE 

RELOCATION 






500 

AVDEN 

REAL 


0 

DATE 

REAL 

ARRAY 

DO 

0 

GP 

REAL 

ARRAY F.P, 

471 

I 

INTEGER 



472 

J 

INTEGER 


0 

K 

INTEGER 

ARRAY 

F.P. 

475 

KC 

INTEGER 


0 

KCOUNT 

INTEGER 

ARRAY 

F.P. 

477 

KCl 

INTEGER 


473 

L 

INTEGER 



465 

M 

INTEGER 


466 

Ml 

INTEGER 



467 

M2 

INTEGER 


470 

M3 

INTEGER 



474 

N 

INTEGER 


0 

NLAT 

INTEGER 


F.P. 

0 

NLONG 

INTEGER 

F.P. 

476 

OZDEN 

REAL 



0 

S 

REAL 

ARRAY F.P. 

0 

SS 

REAL 

ARRAY 

F.P. 



SUBROUTINE ; 

HISTPLT 74/74 OPT-1 FTN 4.7+485 

80/01/23. 

14.08.48, 

1 


SUBROUTINE HISTPLT (ILAT#KK»NLAT#NLONG#RLAT#RKK>RLATBND»RLNGBK> 

RK,HISTPLT 

1 



1NAME,K,L) 

HISTPLT 

2 



COMNON/DD/DATE 

HISTPLT 

3 



DIMENSION ILAT(NLAT)>RLAT(K)»KKCNLAT#NL0N6)#RKK<NLAT#NL0NG) 

HISTPLT 

4 

5 


DIMENSION RLATBNO(K),RLNGBK(L)»RK(L)»NAME(NLAT) 

HISTPLT 

5 



DIMENSION DATE(2) 

HISTPLT 

6 



Kl-K-1 S Ll-L-1 

HISTPLT 

7 



ICOUNT-0 

HISTPLT 

B 



M— 17 

HISTPLT 

9 

10 


NLATl-0 

HISTPLT 

10 



RLATBND(Kl)-0. 

HISTPLT 

11 



RLATBND(K)-1. 

HISTPLT 

12 



RLAT(K1)»0. 

HISTPLT 

13 



RLAT(K)«1. 

HISTPLT 

14 

15 


RLNGBK{L1)«0. 

HISTPLT 

15 



RLNGBK(L)«2. 

HISTPIT 

16 



RK(Ll)-0. 

HISTPLT 

17 



RK(L)-0.1 

HISTPLT 

18 

C 


ILAT CONTAINS THE NUMBER OF DATA POINTS/LATITUDE BAND 

HISTPLT 

19 

20 C 


KK CONTAINS THE NUMBER OF DATA POINTS /GRID BLOCK 

HISTPLT 

20 

C 


FIND MAX VALUES OF ILAT AND KK, ILATMAX AND KKMAX, RESPECTIVELY 

HISTPLT 

21 



ILATMAX-0 $ KKMAX-0 

HISTPLT 

22 

O 


DO 15 I-1,NLAT 

HISTPLT 

23 



IF (ILAT(I).GT. ILATMAX) ILATMAX-ILAT ( I ) 

HISTPLT 

24 

25 


DO 15 J>1,NL0NG 

HISTPLT 

25 



IF (KK(I,J).GT. KKMAX) KKMAX-KK ( I, J ) 

HISTPLT 

26 


15 

CONTINUE 

HISTPLT 

27 



PRINT 101 

HISTPLT 

28 



PRINT 100, ILAT 

HISTPLT 

29 

30 


PRINT 103, DATE 

HISTPLT 

30 



PRINT 104 

HISTPLT 

31 



PRINT 106, DATE 

HISTPLT 

32 



IF (NLAT.LE.18) GO TO 21 

HISTPLT 

33 


19 

ICOUNT-ICOUNT+1 

HISTPLT 

34 

35 


M-M + 18 

HISTPLT 

35 



NLATl-NLATl+lB 

HISTPLT 

36 



PRINT 105, (I,I-M,NLAT1) 

HISTPLT 

37 



DO 20 J»1,NL0NG 

HISTPLT 

38 


20 

PRINT 110, J,(KK(I,J),I«H,NLAT1) 

HISTPLT 

39 

40 


IF (NLAT-IC0UNT+18.GT.18) GO TO 19 

HISTPLT 

40 



M«M+18 

' HISTPLT 

41 


21 

IF (M.LT.O) M-1 

, 1 . ^ ‘ 

HISTPLT 

* . 

42 



4 


I 


SUBROUTINE HISTPLT 74/74 OPT-1 

FTN 

4.7+485 

80/01/23. 

14.06.48 



PRINT 105, (I,I-M,NLAT) 



HISTPLT 

43 



DO 22 J-1,NL0NG 



HISTPLT 

44 

45 

22 

PRINT 110, J,IKK(I,J),I«H,NLAT) 



HISTPLT 

45 


C 

FIND NORNALIZED VALUES OF ILAT AND KK 



HISTPLT 

46 



DO 25 I»1,NLAT 



HISTPLT 

47 


C 

RLAT IS NORMALIZED VALUE OF ILAT 



HISTPLT 

48 



RLAT(I)-ILAT(I) /FLOAT(ILATMAX) 



HISTPLT 

49 

50 


DO 25 J»1,NL0NG 



HISTPLT 

50 


C 

RKK IS NORMALIZED VALUE OF KK 



HISTPLT 

51 



RKK(I, J)-KK(I, J)/FLOAT(KKMAX) 



HISTPLT 

52 


25 

CONTINUE 



HISTPLT 

53 



PRINT 125, ILATMAX,KKMAX 



HISTPLT 

54 

55 


PRINT 126 



HISTPLT 

55 



DO 26 I-1,NLAT 



HISTPLT 

56 


26 

PRINT 127, RLAT(I) 



HISTPLT 

57 



PRINT 103, DATE 



HISTPLT 

58 



PRINT 128 



HISTPLT 

59 

60 


PRINT 106, DATE 



HISTPLT 

60 



ICOUNT-0 



HISTPLT 

61 



M — 17 



HISTPLT 

62 



NLATl-0 



HISTPLT 

63 



IF {NLAT.LE.18) GO TO 31 



HISTPLT 

64 

65 

29 

ICOUNT-ICOUNT+1 



HISTPLT 

65 



M-M+18 



HISTPLT 

66 



NLATl-NLATl+18 



HISTPLT 

67 



PRINT 105, (I,I-M,NLAT1) 



HISTPLT 

68 



DO 30 J«1,NL0NG 



HISTPLT 

69 

70 

30 

PRINT 135, J,(RKK(1, J),I«M,NLAT1) 



HLSTPLI 

TO 



IF (NLAT-IC0UNT*18.GT. 18 ) GO TO 29 



HISTPLT 

71 



M-M+18 



HISTPLT 

72 


31 

IF (M.LT.O) M«1 



HISTPLT 

73 



PRINT 105, (I,I-M,NLAT) 



HISTPLT 

74 

75 


00 32 J«1,NL0NG 



HISTPLT 

75 


32 

PRINT 135, J, <RKK(I, J ),I-M,NLAT) 



HISTPLT 

76 



XL-9. 



HISTPLT 

77 



YL-5. 



HISTPLT 

78 


C 

FOR PROPER SCALING MULTIPLY RLATU) BY 

THE LENGTH 

OF THE Y-AXIS 

HISTPLT 

79 

80 


DO 35 I»1,NLAT 



HISTPLT 

80 



RLATd )-RLAT(I )^YL 



HISTPLT 

81 


35 

RLATBND(I)-(XL/NLAT)/2.+(I-l)*XL/NLAT 



HISTPLT 

82 



DO 36 J-1,NLQNG 



HISTPLT 

83 


36 

RLNGBK(J)>J 



HISTPLT 

84 



SUBROUTINE HISTPLT 7A/7A OPT-1 FTN 'i. 7 + 485 

80/01/23. 

14.06.48 

85 

C 

DX AND DY ARE X AND Y AXES SCALE FACTORS 

HISTPLT 

85 



DX-180./XL 

HISTPLT 

86 



DY-l./YL 

HISTPLT 

87 


C 

T(X OR Y) ■ FREQUENCY OF TIC HARKS/SCALE FACTOR# WHERE SCALE 

HISTPLT 

SB 


C 

FACTOR (DX OP DY) IS A CONVERSION FACTOR BETWEEN AXES UNITS 

HISTPLT 

89 

90 

c 

AND REAL LENGTH 

HISTPLT 

90 



TX— lO./DX 

HISTPLT 

91 



TY — .l/DY 

HISTPLT 

92 


c 

THE INPUT PARAMETERS REQUIRED FOR SUBROUTINE BARPLT HAVE BEEN 

HISTPLT 

93 


c 

CALCULATED. 

HISTPLT 

94 

95 

c 


HISTPLT 

95 


c 

PLOT HISTOGRAMS 

HISTPLT 

96 


c 


HISTPLT 

97 


c 

DATA DISTRIBUTION PER LATITUDE ZONE 

HISTPLT 

98 



CALL NFRAME 

HISTPLT 

99 

100 


CALL AXES(0.»0.# 0 . # X L#-90 . # DX , TX# 0 . # 14HL ATITUDE (DEG)#. 10#- 

14)HISTPLT 

100 



CALL AXES(0.#0.,90.#YL#0.»DY#TY»0.#32HN0RMALIZED NUMBER OF DATA 

POHISTPLT 

101 



1INTS#.10#32) 

HISTPLT 

102 


c 

WBAR IS THE BAR WIDTH WHICH IS ISIZE(SEE DZSTAT PARAMETER LIST) 

HISTPLT 

103 


c 

DEGREES WIDE. 

HISTPLT 

104 

105 


WBAR»(180./NLAT)/DX 

HISTPLT 

105 



CALL BARPLT (RL ATBND# RLAT# NLAT# 1# 1# WBAR#0 ) 

HISTPLT 

106 



CALL NOTATE (2.#6.00# .15»30HNUMBER OF DATA POINTS/LAT BAND#0.#30) HISTPLT 

107 

ro 


CALL NOTATE (2 . > 5. 75# .15# 32HHIST0GRAM INCLUDES DATA FOR DAYS#0. 

#32HISTPLT 

108 



1) 

HISTPLT 

109 

110 


CALL NOTATE ( 2. # 5 . 50# .15# DATE# 0. # 14 ) 

HISTPLT 

110 



ISELECT «0 

HISTPLT 

111 



IF (ISELECT. EQ.O) GO TO 90 

HISTPLT 

IIZ 


c 


HISTPLT 

113 


c 

DATA DISTRIBUTION PER GRID BLOCK 

HISTPLT 

114 

115 


DO 50 I«1#NLAT 

HISTPLT 

115 



DO 45 J-1#NL0NG 

HISTPLT 

116 



45 RK(J)-RKK(l#NLONG+l-J) 

HISTPLT 

117 


c 


HISTPLT 

118 


c 


HISTPLT 

119 

1?0 

c 

NO MODIFICATIONS FOR VARIABLE BLOCK SIZE BELOW THIS POINT. 

HISTPLT 

120 


c 


HISTPLT 

121 


c 


HISTPLT 

122 



CALL NFRAME 

HISTPLT 

123 



CALL AXES (0.#0.#0.#ie.#0.#2.#TX#0.#33HL0NGITUDE INDICES FOR GRID HISTPLT 

124 

125 


lBL0CKS#.15#-33) 

HISTPLT 

125 



CALL AXES(0.#0.#90,#YL#0.#DY#TY#0.#32HN0RMALIZED NUMBER OF DATA 

POHISTPLT 

126 



SUBROUTINE HISTPLT 


74/74 


OPT-1 


FTN 4.7+485 


80/01/23. 14.08.48 



1INTS#.15#32) 

HISTPLT 

127 



CALL NOTATE ( 4 12.> . 1 5, 32HNUMBER OF DATA POINTS/GRID BLOCKS 0.> 32HISTPLT 

128 


1) 

HISTPLT 

129 

130 


CALL NOTATE ( 4. > 11. 75# . 15# 14HLATITUDE INDEX#0.#14) 

HISTPLT 

130 



CALL NOTATE ( 6., 11.75# .15#NAME( I ) #0.# 10) 

HISTPLT 

131 



CALL NOTATE ( 4. # 11. 5# . 15, 32HHIST0GRAM INCLUDES DATA FOR DAYS#0.# 

32HISTPLT 

132 



1) 

HISTPLT 

133 



CALL NOTATE (4 . , 11. 25, . 1 5# DATE# 0 .# 14 ) 

HISTPLT 

134 

135 


CALL BARPLT (RLNGBK#RK# 36# 1# 1#0.25#0) 

HISTPLT 

135 


50 

CONTINUE 

HISTPLT 

136 


90 

CONTINUE 

HISTPLT 

137 



RETURN 

HISTPLT 

138 


100 

format C1X#*ILAT» +#I5) 

HISTPLT 

139 

140 

101 

FORMAT (♦!♦# *THE NUMBER OF DATA POINTS/LATITUDE BAND FROM SOUTH 

THISTPLT 

140 



10 NORTH ARE*) 

HISTPLT 

141 


103 

FORMAT (//*0*#+F0R THE TIME INTERVAL *#A10#A4) 

HlSTiXT 

J.42 


104 

FORMAT (*1*,*NUM8ER OF DATA POINTS/GRID BLOCK*) 

HISTPLT 

143 


105 

FORMAT (*0*#8X#18(I2#5X)/) 

HISTPLT 

144 

145 

106 

FORMAT (1X,A10,A4) 

HISTPLT 

145 


110 

FORMAT (1X#I2#4X#18(I4#3X) ) 

HISTPLT 

146 


125 

FORMAT (*0*#*LATMAX- ♦# 15# 10X#*KKMAX« *#I5) 

HISTPLT 

147 

-i 

126 

FORMAT (*1*, *NORMALIZED NUMBER OF DATA POINTS /LATITUDE BAND*! 

HISIPIT 

148 . 

4^ 

CaI 

127 

FORMAT C1X#*RLAT- *#F7,4) 

HISTPLT 

149 

150 

128 

FORMAT (*1*,*N0RMALIZED NUMBER OF DATA PDINTS/GRID BLOCK*) 

HISTPLT 

150 


135 

FORMAT (1X#I2#4X#18(F5.2#2X) ) 

HISTPLT 

151 



END 

HISTPLT 

152 


SYMBOLIC REFERENCE MAP (R-1) 

FRY POINTS 
3 HISTPLT 

JIABLES SN TYPE RELOCATION 


0 

DATE 

REAL 

ARRAY 

DD 

1140 

DX 

REAL 


L41 

DY 

REAL 



1134 

I 

INTEGER 


L27 

ICOUNT 

INTEGER 



0 

ILAT 

INTEGER 

ARRAY F.P. 

L32 

ILATMAX 

INTEGER 



1145 

ISELECT 

INTEGER 


135 

J 

INTEGER 



0 

K 

INTEGER 

F.P. 



SUBROUTINE AVARPLT 


7A/7^ OPT-1 


FTN ^.7+-!»85 


80/01/23 


14.08.46 


1 


SUBROUTINE AVARPLT < V> A> U#T#M>RAT#STDEVP> STDEVNiNLAT ) 

AVARPLT 

1 



CDMHON/DD/0 

AVARPLT 

2 



DIMENSION V(NLAT)>A(NLAT)#RAT(NLAT)>STDEVP(NLAT)#STDEVN(NLAT) 

AVARPLT 

3 



DIMENSION U(M)»T(M) 

AVARPLT 

4 

5 


DIMENSION D(2) 

AVARPLT 

5 


C 

V AND A ARE THE VARIANCE AND AVERAGE# RESPECTIVELY# OF THE OZONE 

AVARPLT 

6 


C 

DENSITY/LATITUDE BAND 

AVARPLT 

7 


C 

U AND T ARE THE LATITUDE AND OZONE DENSITY INFORMATION/DATA R ECORDAVARPLT 

8 


C 

USED TO PLOT THE SCATTER DIAGRAM. 

AVARPLT 

9 

10 

C 

M IS THE DIMENSION OF U AND T. 

AVARPLT 

10 


c 

STDFVP IS THE STANDARD DEVIATION + MEAN# PROPERLY SCALED TO PLOT 

AVARPLT 

11 


c 

STDEVN IS THE STANDARD DEVIATION - MEAN# PROPERLY SCALED TO PLOT 

AVARPLT 

12 


c 

SUBROUTINE CALPLT (X#Y#IPEN) IS LOCATED IN SECTION 1.4.3 OF 

AVARPLT 

13 


c 

THE GRAPHICS MANUAL 

AVARPLT 

14 

15 

c 

TPEN-2 PEN DOWN 

AVARPLT 

15 


c 

IPEN-3 PEN UP 

AVARPLT 

16 


c 

IPEN LESS THAN ZERO WILL ASSIGN X-O# Y-0 AS THE LOCATION OF 

AVARPLT 

17 


c 

THE PEN AFTER MOVING THE X#Y (CREATE A NEW REFERENCE POINT). 

AVARPLT 

18 


c 

SUBROUTINE PNTPLT ( X# Y# I S YM# IS ) CAN BE FOUND IN SECTION 1.4.70 

AVARPLT 

19 

20 

c 

OF THE GRAPHICS MANUAL. 

AVARPLT 

20 


c 

SUBROUTINE PSEUDO (FN)# FN ■ FILENAME# CAN BE FOUND IN SECTION 

AVARPLT 

21 


c 

1.4.1 OF THE GRAPHICS MANUAL 

AVARPXT 

22 


c 


AVARPLT 

23 


c 

INITIALIZE PARAMETERS 

AVARPLT 

24 

25 


XL-B. 

AVARPLT 

25 



YL-6.0 

0CT79 

1 



DX-180./XL 

AVARPLT 

27 



YMAX-0.65 

0CT79 

2 



YMIN-0.15 

0CT79 

3 

30 


DY»(YMAX-YMIN)/YL 

0CT79 

4 



TX— lO./DX % TY — .1/DY 

AVARPLT 

29 



XT-ABS (TX) 

AVARPLT 

30 



UMAX-0. 

AVARPLT 

31 



UMIN-0. 

AVARPLT 

32 

35 


TMAX — 1. 

0CT79 

5 



TMIN-1. 

0CT79 

6 


c 


AVARPLT 

33 


c 

CONSTRUCT PLOT LABELS AND AXES. 

AVARPLT 

34 



CALL NFRAME 

AVARPLT 

35 

40 


CALL NOTATE ( 2.# 5.00# .15# 15HSCATTER DIAGRAM# 0, #15) 

AVARPLT 

36 



CALL NOTATE (2. #4.75# .15»36HINCLUDES MEAN AND STANDARD DEVIATION# 

OAVARPLT 

37 



1 .#36) 

AVARPLT 

38 



« > I 9 


SUBROUTINE 

AVARPLT 7A/7^ OPT-1. 



FTN A. 

7+A85 

80/01/2 3. 

14.08.^6 



CALL NOTATE ( 2 . 00> A. 50/ . 15# ^ZHDATA 

TAKEN 

FROM 

NIMBUS 

IV BUV MEASURAVARPLT 

39 



1EMENTS»0.,A2) 





AVARPLT 

40 

A5 


CALL NOTATE (2. 00# A. 25/ . 15/ 35HTHIS 

DIAGRAM INCLUDES 

DATA FOR DAYS 

/AVARPLT 

41 



10. #35) 





AVARPLT 

42 . 



CALL NOTATE { 2,00/ A . 00/ . 15/ D/-0. / lA) 





AVARPLT 

43 



CALL AXES (0./0./0./XL/-90./DX/TX/0 

./lAHLATITUDE (DEG )/ , 10/-1A) 

AVARPLT 

44 



CALL AXES (0./0./90./YL/ ,15/DY/TY/O 

./20HT0TAL 

□ZONE 

(ATM-CM)/.10# 

□CT79 

7 

50 


120) 





0CT79 

8 



CALL CALPLT (0./YL/3) 





AVARPLT 

47 



CALL CALPLT (XL/YL/2) 





AVARPLT 

48 - 



CALL CALPLT (XL/0. #2) 





AVARPLT 

49 



CALL CALPLT (0./0./3) 





AVARPLT 

50 

55 


RETURN 





AVARPLT 

51 


C 

« « * * « « * * 1 4> 4c 4c « 4c « % 4c * 4t * 



AVARPLT 

52 



ENTRY SCAT 





AVARPLT 

53 


C 

PLOT SCATTER DIAGRAM - ALSO FIND MAXIMUM 

AND 

MINIMUM 

LATITUDES 

AVARRLT 

54 



DO 30 I-l/H 





AVARPLT 

55 

60 

c 

FIND MAXIMUM AND MINIMUM LATITUDE VALUES 

4c 4c # 4c 4c 4c 4c 4c 4c 4c # 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c A V A R P L T 

56 



IF (Ud) .GT.UMAX) UMAX-U(I) 





AVARPLT 

57 



IF (U(I) .LT.UMIN) UMIN-Ud) 





AVARPLT 

58 


COMMENT — FIND MAXIMUM AND MINIMUM OZONE 

VALUES ^4>4c4c#4c4c«4c4c4cc44c4c4c«4>4c4c4c4c4c 

0CT79 

9 



IF (T(I).GT.TMAX) TMAX»T(I) 





OCT 7.9 

LO 

65 


IF (T(I).LT.TMIN) TMIN-T(I) 





□CT79 

11 


C 

4c 4c 4c 4c 4c 4> 4c 4c 4c 4> 4c 4c 4> >4 4> 4c 4> 4c 4c 4c 4c cf 4c 4c 4> 4c 4> 4< 4c 4c >4 4c 4c % 4c 4c 4c 4c 4c 4c ^ 4>4c 



AVARPLT 

59 


C 

PLOT SCATTER DIAGRAM 





AVARPLT 

60 



X»(U(I)+90.)/DX 





AVARPLT 

61 



Y»(T(I)-YMIN)/DY 





OCT 79 

12 

70 


30 CALL PNTPLT(X/Y/-21/1) 





AVARPLT 

63 



RETURN 





AVARPLT 

64 


C 

% « ^ 4c ^ 4c « 4c 4c 4> 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4 m<c 4c % 4> ♦ 4c 4> 4c « 4c 4c 4c 4c c4 4c 



AVARPLT 

65 



ENTRY ASTO 





AVARPLT 

66 



XS-XL/NLAT 





AVARPLT 

67 

75 


NLATl-NLAT/2 





AVARPLT 

66 



NLAT2-NLAT1+1 





AVARPLT 

69 



DO 35 I-l/NLAT 





AVARPLT 

70 



STDEVP (I )-(A(I )+SORT(V(I ) ) )/DY 





AVARPLT 

71 



STDEVN(I)-(A(I)-SORT(V(I )) )/DY 





0CT79 

13 

80 


STDEVP(I)-STDEVP(I)-YMIN/DY 





0CT79 

14 



STDEVNd )-STDEVN(I )-YMIN/DY 





OCT 79 

15 



35 CONTINUE 





0CT79 

16 



DO AO I-l/NLATl 





AVARPLT 

73 



J-I+NLATl 





AVARPLT 

74 



SUBROUTINE AVARPLT 


7A/74 OPT-1 


FTN A.7+A85 


80/01/23. 14.08.48 


85 

C 


RAT IS THE SCALED DATA MATRIX FOR SPACING PLOTTED AVERAGES AND 

AVARPLT 

75 


C 


STANDARD DEVIATION ALONG THE X - AXIS 

AVARPLT 

76 




RAT(I>»XS/2.+(NLAT/2-l+I)+XS 

AVARPLT 

77 




RAT(J)«XS/2.+(NLAT-J)*XS 

AVARPLT 

78 



40 

CONTINUE 

AVARPLT 

79 

90 

COMMENT — FIND LATITUDE INDEXES# IMAX AND IMIN# CORRESPONDING TO 

AVARPLT 

80 


C 


UMAX AND UMIN. THEN# CALCULATE AN ADJUSTED VALUE OF • 

AVARPLT 

81 


C 


RAT(IMAX) AND RAT(IHIN) SUCH THAT THE EXTREME MEANS WILL 

AVARPLT 

82 


C 


BE PLOTTED IN THE END LATITUDE ZONES HALF-WAY BETWEEN THE 

AVARPLT 

83 


C 


ZONE"S BEGINNING AND THE EXTREMUM LATITUDE VALUES. 

AVARPLT 

84 

95 



IMAX-UHAX/(180/NLAT) 

AVARPLT 

85 




LMAX-IMAX*(180/NLAT) 

AVARPLT 

86 




IMAX-IMAX+1 

AVARPLT 

87 


COMMENT — RATMAX IS THE HALF-WAY POINT FDR THE EXTREME MAXIMUM 

AVARPLT 

88 


C 


LATITUDE ZONE. 

AVARPLT 

89 

100 



RATMAX-(UMAX-LMAX)/ (2.+DX) 

AVARPLT 

90 




RAT { IMAX )»RAT( IMAX >-XS/2. 

AVARPLT 

91 




RAT (IMAX )-RAT(IMAX)+RATMAX 

AVARPLT 

92 




IMIN-UMIN/(180/NLAT) 

AVARPLT 

93 




LMIN-IMIN*(180/NLAT) 

AVARPLT 

94 

__105 



IMIN»(NLAT/2+l)-IMIN 

AVARPLT 

95 

4::> 

fT\ 

COMMENT — RATMIN IS THE HALF-WAY POINT FOR THE EXTREME MINIMUM 

AVARPLT 

96 

U 1 
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LATITUDE ZONE. 

AVARPLT 
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RATMIN-(UMIN-LMIN)/(2.*DX) 

AVARPLT 

98 




RAT(IMIN)«RAT{IMIN)+XS/2. 
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99 

110 



RAT{IMIN)-RATaMIN)+RATMIN 
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4c 4. « ^ 4t ^ 4: % 4c * 4: * * « * « 
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DO 41 I-1#NLAT 

AVARPLT 
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A NOW BECOMES THE PROPERLY SCALED AVERAGE TO BE PLOTTED 

AVARPLT 
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ALONG THE Y-AXIS 
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41 

A(I)-(A(I)-YMIN)/DY 
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4 4 1 # * 4< ^ 4< ♦ * 4> « 1# 4c 4c 4< 4> ^ <4 4 m^ 4c 4i # ♦ 4 4c 4c 4> 4i 4> 4c 4c 4> 4c 4c 4c 4c 4> « 4c 4< 
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NOW PLOT MEANS AND DRAW IN CONNECTING "CURVE” 
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IF (A(1).E0.-YMIN/DY) GO TO 42 
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CALL PNTPLT (RAT( D# A( 1 )»-ll#2 ) 

AVARPLT 

108 

120 


42 

CONTINUE 

DCT79 
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DO 45 I-2#IMAX 

AVARPLT 
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IF (Ad-D.EO.-YMIN/DY) GO TO 43 

0CT79 
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CALL CALPLT (RAT( I)#A (I) #2) 
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110 



43 

CONTINUE 

0CT79 
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IF (A(I).EQ.-YMIN/DY) GO TO 45 
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CALL PNTPLT { RATI I ) # A ( I ) #-11# 2 ) 
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SUBROUTINE AVARPLT 


74/7A OPT-1 


FTN A. 7+485 


80/01/23. 


14.08.48 
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135 


140 


145 


150 


155 


160 


165 


C 


C 


C 


45 

CONTINUE 




0CT79 

24 


CALL CALPLT ( RAT( 1 ) > A (1 ) >3 ) 
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DO 50 I«NLAT2#IMIN 




AVARPLT 
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IF (A(l) .EQ.-YHIN/0Y.AND.I.E0.NLAT2) GO 

TO 

46 


OCT 79 

25 


IF (A(I).EQ.-YMIN/DY.0R.A(I-1».EQ.-YHIN/DY. 

AND. I 

.GT.NLAT2) 

GOTO 460CT79 

26 


CALL CALPLT (R AT ( I ) # A ( I ) , 2 ) 




AVARPLT 
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46 

CONTINUE 




0CT79 

27 


IF (Ad) .EQ.-YMIN/DY) GO TO 50 




DCT79 

28 


CALL PNTPLT ( RAT ( I ) > A (I ) *-ll»2 ) 




0CT79 

29 

50 

CONTINUE 




OCT 79 

30 


PLOT STANDARD DEVIATION 
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DO 55 I«1#ININ 




AVARPAT 

117 


IF {(I.GT.IMAX).AND.(I.LT.NLAT2)) GO TO 

55 
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IF (Ad). EQ.-YMIN/DY) GO TO 55 
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CALL CALPLT ( ( RAT ( I )-0. 06 ) ,S TDE VN (I )^ 3 ) 




AVARPLT 

119 


CALL CALPLT ( ( RAT ( 1 ) +0. 06)# STDEVN (1)> 2 ) 
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CALL CALPLT (RAT (I )# STDE VN (I )# 3 ) 
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CALL CALPLT (RATd )#STDEVP(I )>2) 
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CALL CALPLT ( ( RAT (I )-0. 06) # STDEVP (I )# 3 ) 
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123 


CALL CALPLT ( ( R AT (I ) + 0.06 ) ,STDE VP ( I )# 2 ) 




AVARPAT 

124 

55 

CONTINUE 
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125 


PLOT MEANS AND STANDARD DEVIATIONS AS A 

SEPERATE 

FRAME. 

AVARPAT 

1.26 


CALL NFRAME 
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127 

CALL NOTATE ( 2 .00# 5. 00# . 15# 27HMEAN AND STANDARD DEVIATI0N#0.#27) 
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126 

CALL NOTATE (2. 00#4.75, . 15# 42HDATA TAKEN FROM 

NIMBUS 

IV BUV MEASURAVARPIT 

129 

1EMENTS#0.#42 ) 



AVARPLT 

130 

CALL NOTATE (2 . 00# 4. 5 0, . 15# 35HTHIS DIAGRAM INCLUDES 

DATA FOR DAYS#AVARPLT 

131 

10. #35) 



AVJlRPLT 

032 

CALL NOTATE (2 . 00# 4. 25# . 1 5# D#0 .# 14 ) 
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133 

CALL AXES (0.#0.#0.#XL#-90.#DX#TX#0.#14HLATITUDE ( DE G ) # . 10#-14 ) 

AVARPLT 

134 

CALL AXES (0.#0.#90.#YL#.15#DY#TY#0.#20HTOTAL 

OZONE 

(ATM-CM)#.10# 

QCT79 

32 

120) 
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33 

CALL CALPLT (0.#YL#3) 



AVARPLT 

137 

CALL CALPLT (XL#YL#2) 



AVARPLT 

138 

CALL CALPLT (XL#0.,2) 
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139 

CALL CALPLT (0.#0.#3) 



AVARPLT 

140 

NOW PLOT MEANS AND DRAW IN CONNECTING ’’CURVE” 



AVARPLT 

141 

IF (A(l) .EQ.-YMIN/DY) GO TO 556 
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CALL PNTPLT (RAT( 1 )#A (1 )#-22#2 ) 
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35 

556 CONTINUE 
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36 

DO 56 I»2#IMAX 



AVARPLT 
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IF (Ad-l). EQ.-YMIN/DY) GO TO 557 



0CT79 

37 



SUBROUTINE 

AVARPLT 7A/7A OPT-1 FTN A.7+A85 

80/01/23. 

lA.Oe.Af 



CALL CALPLT (RAT ( I ) > A ( I J #2 ) 

AVARPLT 

lAA 

170 

557 

CONTINUE 

0CT79 
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IF (A(I) .EQ.-YMIN/DY) GO TO 56 

0CT79 

39 



CALL PNTPLT (RAT(I )#A(I )#-22#2) 

0CT79 

AO 


56 

CONTINUE 

0CT79 
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CALL CALPLT (R AT ( 1 )« A ( 1 ) » 3 ) 

AVARPLT 

1A6 

175 


DO 57 I-NLAT2>ININ 

AVARPLT 

1A7 



IF (A(l) .EQ,-YMIN/DY,AND.I.EQ.NLAT2) GO TO 558 

0CT79 

A2 



IF CA(I) .E0.-YMIN/DY.0R.A(I-1),EQ.-YMIN/DY.AND.I.GT.NLAT2) G0T0558DCT79 

A3 



CALL CALPLT (RAT( I) j A( I ) » 2 ) 
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1A8 


558 

CONTINUE 

0CT79 
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IF (A(I) .EQ.-YMIN/DY) GO TO 57 
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CALL PNTPLT (R AT (I ) # A ( I ) #-22# 2 ) 

0CT79 

A6 


57 

CONTINUE 

0CT79 
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PLOT STANDARD DEVIATION 
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DO 58 I-1#IMIN 

AVARPLT 

151 

185 


IF (d.GT.IMAX).AND.(I.LT.NLAT2)) GO TO 58 
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IF (A(I). EQ.-YMIN/DY) GO TO 58 

0CT79 
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CALL CALPLT ( ( RAT (I )-0. 06) #STDEVN (I ), 3 ) 

AVARPLT 

153 



CALL CALPLT ( ( R AT ( I ) +0.06)# STDE VN( I ) » 2 ) 

AVARPLT 

15A 



CALL CALPLT ( R AT( I ) # STDE VN ( I )# 3 ) 

AVARPLT 

155 

190 


CALL CALPLT ( R AT ( I ) #STDEVP ( I )# 2 ) 

AVARPLT 
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CALL CALPLT ( ( RAT ( I )-0. 06) # STDE VP( I )# 3 ) 

AVARPLT 

157 



CALL CALPLT ( ( R AT (I )+0. 06 ) #STDE VP ( I )# 2 ) 

AVARPLT 
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58 

CONTINUE 
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PRINT 110 

AVARPLT 

160 

195 


PRINT 115# D 

AVARPLT 

161 



DO 60 I-1#NLAT 

AVARPLT 

162 


60 

PRINT 100, I,PAT(I),A(I),V(I),STDEVP(I)#STDEVN(I) 

AVARPLT 

163 



PRINT 105# M,UMAX#UMIN 
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PRINT 120#TMIN,TMAX 

0CT79 

A9 

200 

120 

FORMAT (1X,*TMIN« ♦ , E15 . 8, 5X,*TMAX» *,E15.8) 

0CT79 
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RETURN 

AVARPLT 
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100 

FORMAT (♦0*#*I« *#I2#<tX, ♦RAT* ♦# F6. 2, <»X# *A» ♦# F6.2# AX#*V* *,£11. 

A, AVARPLT 
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1AX#*ST0EVP» ♦#F6.2,<rX#*STDEVN« *#F6.2) 
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105 

FORMAT (♦0^,^M- ♦,I6,5X,*UMAX- ♦,F7.3,5X,*UMIN- *,F7.3) 
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110 

FORMAT (♦!♦, ♦X-AXIS SCALE -RAT-, AVERAGES, VARIANCES# AND STANDARD 

DEAVARPLT 

169 


IVIATIONS USED IN AVARPLT*) 

AVARPLT 

170 


115 

FORMAT (T9, ♦FOR THE TIME PERIOD ♦,A10,A^) 
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END 
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